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Introduction 


Described  herein  are  the  key  findings  and  results  associated  with  the  project  entitled 
“Investigation  of  ELF  Signals  Associated  with  Mine  Warfare,  A  University  of  Idaho  and 
Acoustic  Research  Detachment  Collaboration,  Phase  Two.”  Phase  Two  is  a  continuation 
of  the  Phase  One  effort  under  the  same  title.  The  scope,  objectives  and  outcomes  of  Phase 
Two  are  similar  to  those  described  in  the  report  and  proposal  of  Phase  One.  Some  of  the 
following  text  is  also  found  in  the  Phase  One  report. 

Extremely  low  frequency  (ELF)  electromagnetic  signals  are  used  by  enemy  combatants 
to  detect  and,  subsequently,  to  incapacitate,  by  means  of  surface  and  subsurface  mines, 
naval  vessels.  This  program  is  of  high  importance  to  the  Navy  -  particularly  since  ELF 
signals  are  one  of  the  primary  signature  emissions  of  the  Navy's  proposed  electric  ship 
fleet. 

The  questions  that  are  being  asked  in  this  investigation  are:  1)  once  an  ELF  signal  is 
generated,  how  far  will  it  propagate  and  still  be  detectable  and  2)  how  can  such  signals  be 
modeled,  excited  and  measured?  To  this  end,  the  scenario  considered  is  one  in  which  an 
ELF  source  of  the  electric  or  magnetic  kind  is  located  in  or  above  water,  such  as  a  lake  or 
ocean.  This  source  stimulates  an  ELF  signal  that  is  free  to  propagate  in  the  water  and  air, 
and  is  reflected  by  various  material  interfaces,  say  between  the  water  and  air,  or  between 
the  water  and  the  floor.  For  purposes  of  experimental  demonstration,  the  investigation 
focuses  on  the  scenario  of  ELF  sources  and  signals  in  the  context  of  Lake  Pend  Oreille, 
where  the  Acoustic  Research  Detachment  (ARD,  Bayview,  Idaho)  is  located  and 
entrusted  with  the  necessary  assets  to  perform  validation  measurements. 

The  research  program  was  designed  with  two  major  thrusts:  Modeling  and 
experimentation.  The  modeling  thrust  was  coordinated  and  executed  by  the  University  of 
Idaho  (UI),  Moscow,  Idaho;  the  experimentation  thrust  was  coordinated  and  executed  by 
ARD.  This  report  focuses  primarily  on  the  modeling  thrust.  A  separate  report  from  ARD 
has  been  issued  that  addresses  the  experimentation  thrust  (See  “ELF  Phase  Three  Test  1,” 
complied  by  Frank  Jurenka,  Chris  Burgy  and  Vicki  Pfeifer,  July  15,  2010.  Note:  ARD 
expensed  Phase  Three  funds  while  UI  was  still  expensing  Phase  Two  funds.). 

Both  students  and  faculty  of  the  University  of  Idaho  and  of  Washington  State  University 
were  involved  in  this  project.  Team  members  include: 

Prof.  Jeffrey  L.  Young  (UI),  Lead  PI: 

•  Dr.  Christopher  L.  Wagner,  Research  Engineer,  FDTD  code  development 

•  Mr.  Robert  Rebich,  MSEE  RA,  Quasi-static  code  development 

•  Mr.  Christopher  Johnson,  MSEE  RA,  Data  analysis  and  code  development 
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•  Mr.  Das  Butherus,  MSEE  RA,  Data  analysis  and  code  development 

•  Mr.  Chenchen  “Jimmy”  Li,  BSEE  RA,  Topographical  data  translation 

•  Ms.  Neelima  Dahal,  BSEE  RA,  Data  analysis 

•  Mr.  Markus  Geiger,  BSEE  RA,  Data  analysis 
Prof.  Dennis  Sullivan  (Ul): 

•  Mr.  Yang  Xia,  Research  engineer,  FDTD  parallelization  code  development 

•  Mr.  Alireza  Mansoori,  MSEE  RA,  FDTD  parallelization  code  development 
Prof.  Robert  Olsen  (WSU): 

•  Mr.  Zhi  Li,  MSEE  RA,  Layered  media  modeling 

ELF  Modeling 

The  activities  pursued  during  Phase  Two  continue  those  pursued  during  Phase  One,  with 
particular  emphasis  on  the  refinement  and  validation  of  the  numerical  models.  Portions  of 
the  following  text  summarize  the  modeling  effort  undertaken  by  the  University  of  Idaho 
and  are  also  found  in  the  Phase  One  report.  Some  of  that  text  has  been  updated  to  reflect 
new  knowledge  gained  since  the  Phase  One  report  was  issued.  In  subsequent  sections, 
unique  results  and  findings  associated  with  Phase  Two  activities  will  be  presented. 

Modeling  of  ELF  electromagnetic  signals  in  water  environments  can  be  accomplished 
either  by  means  of  direct,  analytical  solution  of  Maxwell's  equations  or  by  numerical 
solutions  of  the  same.  The  former  is  attractive  for  purposes  of  gaining  insights  into  the 
physical  mechanisms  that  hinder  or  aid  the  propagation  of  ELF  signals.  The  disadvantage 
is  found  in  the  number  of  simplifying  assumptions  that  are  made  to  bring  about  a 
closed-form  solution.  A  numerical  solution  has  no  such  simplifying  assumptions,  but  does 
suffer  from  discretization  errors.  In  principle,  it  can  model  all  of  the  physical  and 
geometrical  features  of  the  domain  of  consideration.  The  price  paid  for  doing  so,  however, 
is  the  required  time  and  the  CPU/memory  resources  needed  to  accomplish  the  task.  Data 
visualization  and  management  are  other  issues  that  need  to  be  addressed  when  working 
with  large  data  sets  produced  by  numerical  solvers.  The  positive  and  negative  tradeoffs 
between  these  two  approaches  (i.e.  analytical  vs.  numerical)  suggest  that  no  one  method 
is  superior.  For  that  reason  the  UI  team  adopted  a  diverse  strategy  that  encompasses  many 
different  approaches  in  order  to  assure  a  positive  outcome  and  to  provide  deliverable 
modeling  methodologies. 

The  five  principle  techniques  or  tools  that  were  considered  during  the  Phase  One  and 
Phase  Two  efforts  were  the  a)  Sommerfeld  Half-Space  (SHS)  method,  b) 
Finite-difference,  time-domain  method  (FDTD),  c)  High  Frequency  Structural  Simulator 
(HFSS),  finite-element  code,  d)  Maxwell  code  and  e)  quasi-static  method  (QES).  A 
summary  of  these  methods  is  provided  next.  Detailed  technical  information  on  the  SHS, 
FDTD  and  quasi-static  methods  are  provided  in  the  attached  appendices. 
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Sommerfeld  Half-Space  (SHS)  Method 


The  SHS  method  is  an  analytical  approach  that  assumes  that  all  interfaces  (say  between 
water  and  air,  or  between  water  and  floor)  are  planar  and  infinitely  extended.  This 
assumption  is  reasonably  valid  for  the  water-air  interface,  particularly  in  open  water 
regions  where  the  source  is  located  near  the  surface.  For  the  littoral  zones,  the  method 
may  fail,  particularly  when  electric  sources  are  used  to  excite  the  ELF  signals.  By 
assuming  that  the  interfaces  are  flat,  a  closed-form  solution  can  be  devised  that  is  cast  in 
terms  of  Fourier-Bessel  integrals.  These  integrals  can  be  evaluated  numerically  and 
rapidly  in  a  matter  of  seconds  on  any  desktop  machine.  Even  with  the  potential 
deficiency  of  treating  all  interfaces  as  planar,  the  SHS  method  is  attractive  as  a  validation 
tool  for  the  other  numerical  modeling  approaches.  For  example,  the  team  used  the  SHS 
method  to  validate  the  data  produced  by  the  FDTD  or  HFSS  methods  (described  next) 
when  these  numerical  methods  consider  the  same  layered  media  problem  statement.  The 
SHS  method  is  also  attractive  in  quantifying  the  up-over-down  effect.  This  effect  is 
associated  with  a  low  signal  loss  path  through  the  air  and  a  high  signal  loss  path  up  and 
down  through  the  water.  If  the  path  through  the  water  is  short,  then  the  up-over-down 
signal  loss  can  be  low  relative  to  a  direct  path  between  a  source  and  sensor  in  the  water. 
Professor  Robert  Olsen  of  Washington  State  University  (WSU)  is  the  lead  investigator  of 
the  SHS  method. 

Finite  Difference.  Time-Domain  Method  (FDTD) 


The  FDTD  method  is  a  numerical  approach  that  discretizes  Maxwell's  equations  in  their 
fundamental  form  using  a  staggered  grid  and  leap-frog  integrator.  This  method  has  been 
fully  vetted  in  the  open  literature  and  has  been  established  as  a  robust  way  of  obtaining 
accurate  simulation  data.  In  principle,  the  FDTD  method  accounts  for  all  material 
interfaces  and  inhomogeneities  by  assigning  permittivity,  permeability  and  conductivity 
values  along  edges  of  the  grid  elements.  Curvilinear  boundaries  are  approximated  by 
straight  line,  stair-stepped  boundaries.  For  geometrical  features  that  are  significantly  less 
than  a  wavelength,  such  stair-stepping  causes  no  appreciable  errors  in  the  computed  data. 
Note  that  the  domain  of  interest  at  Lake  Pend  Oreille  does  not  exceed  8  km  on  a  side;  the 
lake  floor  at  its  deepest  point  is  about  335  m.  Assuming  an  operating  frequency  of  100  Hz 
and  a  water  conductivity  of  0.018  S/m,  we  note  that  the  corresponding  skin  depth  is  375 
m  and  the  wavelength  is  2.356  km;  for  air,  the  wavelength  is  3,000  km.  Thus  the  domain 
spans  a  fraction  of  a  wavelength  in  air  but  about  3.4  wavelengths  (or  21.3  skin  depths)  in 
water.  The  significant  disparity  between  these  two  relative  sizes  potentially  introduces 
computational  complexities.  One  area  of  concern  is  the  proper  design  of  an  absorbing 
boundary  condition  (ABC)  or  perfectly  matched  layer  (PML)  that  will  allow  an  open 
physical  domain  to  be  truncated  into  a  finite  computational  domain.  Placement  of  this 
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ABC/PML  in  terms  of  wavelengths  is  critical  if  non-spurious  reflections  are  to  be 
avoided.  Significant  time  and  effort  was  expended  to  figure  out  a  way  to  design  an 
optimal  PML.  The  outcome  of  this  effort  is  described  in  the  paper  by  Wagner  and  Young, 
“FDTD  numerical  tests  of  the  convolutional-PML  at  extremely  low  frequencies,”  IEEE 
Antennas  and  Wireless  Propagation  Letters,  vol.  8,  pp.  1398-1401,  2009.  This  paper  is 
attached  to  this  report.  Dr.  Christopher  Wagner  of  the  University  of  Idaho  is  the  lead 
researcher  of  the  FDTD  and  PML  effort. 

Given  the  amount  of  computing  time  that  is  required  to  run  an  FDTD  simulation  for 
domain  sizes  contemplated  in  this  project  (domains  are  on  the  order  of  kilometers), 
simulation  times  can  be  excessive  (i.e.  hours  to  days).  For  this  reason,  an  effort  in 
computational  parallelization  was  undertaken  using  message  passing  interface  (MP1) 
protocols  and  specialized  graphics  hardware.  As  shown  in  Appendix  C,  simulation  times 
can  be  reduced  by  over  a  factor  16  using  these  kinds  of  parallelization  techniques. 
Professor  Dennis  Sullivan  of  the  University  of  Idaho  is  the  lead  researcher  of  this  activity. 

High  Frequency  Structural  Simulator  (HFSS) 

HFSS  is  a  commercially  available  electromagnetic,  finite-element,  frequency-domain, 
numerical  solver  that  has  been  designed  by  Ansoft/Ansys  for  antenna  and  microwave 
circuit  applications.  One  question  that  was  asked  in  this  investigation  was  whether  such  a 
tool  could  be  used  to  predict  the  electromagnetic  propagation  characteristics  of  an  ELF 
signal  in  a  highly  conductive  environment.  In  Phase  One  and  Two,  the  answer  to  this 
question  was  inconclusive  due  to  source  modeling  issues.  (In  the  Phase  Three  report, 
however,  the  answer  will  be  more  conclusive  and  positive.)  Professor  Jeffrey  Young  of 
the  University  of  Idaho  is  the  lead  engineer  of  the  HFSS  effort. 

Maxwell 


Maxwell  is  also  a  commercial  code  developed  by  Ansoft/Ansys.  However,  unlike  HFSS, 
it  is  a  static  solver  for  either  electric  or  magnetic  fields.  Since  ELF  waves  are  static-like 
in  the  vicinity  of  the  source,  questions  that  have  been  raised  by  the  team  are  these:  1)  At 
what  distance  are  the  fields  more  static-like  rather  than  wave-like  and  2)  can  ELF  waves 
be  modeled  by  a  static  solver  in  some  region  about  the  source.  Professor  Jeffrey  Young  of 
the  University  of  Idaho  is  the  lead  engineer  of  this  effort. 

Quasi-Static  Method 

A  custom  quasi-static  method  was  also  considered  given  that  ELF  signals  are  quasi-static 
in  the  vicinity  of  the  source.  By  definition,  the  quasi-static  method  does  not  consider  any 
wavelike  mechanisms  in  Maxwell’s  equations;  it  assumes  that  the  field  lines  are  the  same 
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as  the  static  field,  but  oscillating.  This  is  accomplished  by  neglecting  magnetic 
displacement  currents  for  electric  sources.  By  doing  so,  simple  solutions  can  be 
constructed  that  correlate  well  with  other  more  advanced  solutions,  like  HFSS  and  FDTD. 
Professor  Jeffrey  Young  of  the  University  of  Idaho  is  the  lead  engineer  of  the  quasi-static 
modeling  effort. 

Lake  Parameters  and  Discretization 


Unlike  the  December  2008  experiments  in  which  the  experiments  were  conducted  in  an 
open  area  of  the  lake,  the  domain  for  the  March  2010  experiment  encompasses  significant 
geometrical  features  above  and  below  water.  This  was  purposefully  chosen  to  be  so  in 
order  to  exercise  the  limits  of  the  various  numerical  and  analytical  models.  That  is,  we 
would  expect  the  ELF  signals  in  the  open  area  to  be  far  easier  to  model  than  those  in  a 
more  cluttered  environment,  due  to  the  changes  in  the  material  parameters  and 
geometrical  features  of  the  environment.  Hence,  we  wanted  the  most  severe  environment 
possible  to  see  if  the  models  would  fail  to  produce  the  correct  data. 

The  FDTD  and  HFSS  numerical  methods  require  a  precise  understanding  of  the  electrical 
and  geometrical  features  of  the  lake.  The  domain  of  interest  considered  in  Phase  Two  is 
the  area  known  as  ldlewilde  Bay  and  is  shown  below.  The  domain  is  about  6  km  by  7  km 
on  a  side  and  represents  the  general  area  where  actual  experiments  were  performed  in 
March  2010  using  both  electric  and  magnetic  sources.  (See  the  report  entitled  “ELF 
Phase  Three  Test  1”  complied  by  Frank  Jurenka,  Chris  Burgy  and  Vicki  Pfeifer,  July  15, 
2010.  Note:  ARD  expensed  Phase  Three  funds  while  U1  was  still  expensing  Phase  Two 
funds.) 


The  terrain  elevation  data  (relative  to  sea  level)  along  with  their  corresponding 
coordinates  (in  varying  forms)  were  extracted  from  three  sources:  a  data  set  from 
insideidaho.org,  an  AUTOCAD  file  of  Lake  Pend  Oreille  Contours  from  the  Idaho 
Geological  Survey,  and  data  points  taken  manually  from  a  provisional  map  of  Lake  Pend 
Oreille.  The  coordinates  of  each  data  point  were  converted  into  meters  northing  and 
easting  in  Idaho  West  State  Plane;  any  elevation  data  in  feet  were  converted  to  meter  - 
thus,  all  three  data  sets  conform  to  the  same  system.  All  three  data  sets  were  compiled 
together  (minor  adjustments  were  made  to  eliminate  conflict  between  the  data  sets). 
Interpolation  of  elevation  data  at  all  points  along  two  vectors  (in  x  and  y  direction  that 
define  the  area  to  interpolate)  was  accomplished  using  the  'griddata'  function  in  Matlab. 
This  created  a  matrix  height  field  that  defines  the  elevation  and  depth  of  the  terrain  or 
lake  at  each  point  in  lm  intervals.  The  matrix  height  field  was  then  used  as  an  input  file 
for  the  various  numerical  solvers,  i.e.  FDTD,  HFSS  or  Maxwell. 

The  height  field,  if  used  with  HFSS  or  Maxwell,  needs  to  be  converted  into  a  solid  model. 
The  first  step  is  to  extract  data  from  the  height  field  into  x,y,z  coordinates.  Then,  in 
AutoCAD,  the  command  '3dmesh'  is  used  to  create  a  mesh  that  is  up  to  255x255  cells  in 
dimension  from  those  coordinates.  Since  the  height  field  is  6240x7520  cells  in  size,  the 
data  is  down-sampled  so  that  it  will  be  within  the  bounds  of  '3dmesh'.  After  meshing,  an 
AutoCAD  script,  'M2S-2007.1sp'  is  used  to  convert  the  mesh  into  a  solid  figure.  This 
solid  figure  is  then  exported  as  an  ACIS  .sat  file  (which  is  supported  by  HFSS).  However, 
the  mesh  on  the  surface  of  the  solid  is  too  refined  and  uniform  for  HFSS  to  use  efficiently 
in  data  computation.  Therefore,  an  additional  remeshing  step  is  necessary  via  the  mesh 
tool  Cubit.  By  combining  all  the  surfaces  of  the  original  mesh  into  one  composite  surface, 
the  composite  surface  is  then  meshed  using  one  of  Cubit's  meshing  schemes. 
Unfortunately,  Cubit  cannot  imprint  the  new'  mesh  onto  the  original  AutoCAD  solid;  the 
new  mesh  must  be  converted  into  a  solid  itself.  The  mesh  is  exported  into  an  .inp  file  and 
then  re-imported  into  Cubit,  which  removes  the  AutoCAD  solid  and  leaves  only  the  Cubit 
mesh.  The  mesh  is  then  converted  into  a  solid  within  Cubit  and  is  exported  back  into  an 
ACIS  .sat  format.  Clearly,  this  is  an  involved  process,  but  a  necessary  one  when  using 
HFSS  or  Maxwell. 

In  addition  to  precise  geometrical  data,  the  various  solvers  also  require  precise 
knowledge  of  the  conductivity  of  the  lake  and  the  mud  at  the  bottom  of  the  lake.  The  UI 
team  used  a  value  of  0.018  S/m  for  the  water  and  0.012  S/m  for  the  mud  floor.  These 
numbers  were  previously  measured  by  ARD  during  Phase  One.  As  for  the  value  of  the 
dielectric  permittivity  of  the  lake,  this  was  not  deemed  essential,  since  displacement 
currents  in  the  lake  are  virtually  insignificant  relative  to  the  conduction  currents. 

It  should  be  noted  that  a  major  shortcoming  of  the  modeling  effort  has  nothing  to  do  with 
the  modeling  methodology,  but  with  the  lack  of  information  about  the  environment  to  be 
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modeled.  For  example,  we  treat  the  problem  statement  as  if  the  environment  is  only 
comprised  of  three  homogeneous  substances:  water,  air  and  mud.  Clearly,  this  is  not  so. 
The  lake  bottom,  which  we  call  mud,  is  actually  an  inhomogeneous  substance  of  rock  and 
silt  that  is  saturated  by  water.  The  land,  which  is  called  mud,  is  an  inhomogeneous 
substance  of  rock,  dirt,  trees  and  structures.  Only  the  water  and  air  are  homogeneous  for 
which  numbers  like  permittivity  and  conductivity  are  known.  Hence,  errors  between 
experimental  data  and  simulation  data  can  be  attributed  to  the  lack  of  knowledge  of  the 
environment  and  certain  guesses  about  the  quantification  of  the  environment. 

Electric  and  Magnetic  Sources 

Two  kinds  of  electric  sources  were  used  in  the  March  2010  experiment:  1)  a  4  meter,  2 
Ampere  (max)  electric  source  placed  on  a  boat  hull  that  skimmed  the  surface  of  the  water 
and  2)  a  15  meter,  3  Ampere  (max)  portable  electric  source  that  was  lowered  from  15 
meters  in  the  water  to  the  lake  floor  (i.e.  about  152  meters).  Additionally,  a  3.6  meter  by 
3.6  meter,  12  turn,  20  Ampere  magnetic  source  was  also  used  to  stimulate  ELF  signals; 
this  source  was  rigidly  placed  on  the  shore  at  Farragut  State  Park.  For  both  electric  and 
magnetic  sources,  the  ELF  signals  were  measured  using  a  portable  electromagnetic  array 
(EMA)  that  was  lower  into  the  water  at  depths  ranging  from  15  m  to  152  m.  Source  and 
sensor  locations  associated  with  the  March  2010  experiment  are  shown  below. 
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Experimental  Data  Post-Processing 


Once  the  experimental  data  was  collected,  it  was  transferred  to  the  University  of  Idaho  in 
binary  format  for  processing  -  particularly,  to  extract  the  desired  frequency  domain 
signals  from  the  time-domain  data.  The  first  step  was  to  pre-process  the  measured  data 
into  a  useable  format  and  to  scale  the  data  using  appropriate  scaling  factors,  as  provided 
by  ARD.  Next  the  data  was  analyzed  and  plotted  to  identify  experimental  runs  that 
correspond  to  fairly  stationary  source  locations,  since  the  ELF  models  assume  both 
stationary  sources  and  sensors.  Typically,  ten  seconds  of  contiguous  data  sets  can  be 
obtained,  which  are  transformed  into  the  frequency  domain  using  fast  Fourier  Transform 
(FFT)  methods.  The  FFT  data  reveal  the  time-harmonic  signal  strength  of  each 
electromagnetic  field  component  relative  to  the  coordinate  system  of  the  experiment.  To 
compare  these  data  with  simulation  data,  coordinate  transformations  are  performed  on  the 
experimental  data  using  GPS  sensor  data.  The  processing  of  the  data  is  not  completely 
automatic  given  random  noise  spikes  and  discontinuities  in  the  data  streams.  To  avoid 
post-processing  conversion  anomalies,  only  clean,  contiguous  data  were  processed. 

Results 

Data  Comparisons:  Simulation  Methods  vs  Sommerfeld  Half  Space  Method 

As  noted  previously,  the  Sommerfeld  Half-Space  (SHS)  Method,  being  an  exact  solution 
of  Maxwell’s  equations,  can  be  used  to  benchmark  the  accuracy  of  the  various  methods 
employed  in  this  project.  Two  sets  of  plots  are  shown  below  for  vertical  magnetic  dipole 
excitation  and  horizontal  electric  dipole  excitation.  For  the  former,  it  is  clearly  seen  that 
the  Sommerfeld  data  (identified  as  WSU,  who  were  the  developers)  and  the  FDTD  data 
are  closely  correlated,  thus  validating  the  FDTD  methodology  and  code.  There  is  fairly 
good  correlation  between  the  data  sets  of  WSU  and  Maxwell  and  no  correlation  between 
the  data  sets  of  WSU  and  HFSS.  This  poor  correlation  is  attributed  to  the  way  HFSS 
models  Hertzian  dipoles  in  lossy  media.  Although  there  are  some  “tricks”  for  getting 
better  data,  these  tricks  involve  the  use  of  scaling  factors  that  cannot  be  rigorously 
justified  by  theory.  Moreover,  a  priori  reliability  is  never  assured.  (Research  conducted 
during  Phase  Three  has  found  much  more  reliable  ways  to  assure  good  data;  this  will  be 
reported  in  the  Phase  Three  report.)  Similar  conclusions  can  be  reached  for  the  vertical 
electric  dipole  case,  but  with  additional  validation  of  the  quasi-static  method.  Since 
Maxwell  only  predicts  electric  fields  by  assuming  no  excitation  of  the  magnetic  field  and 
since  the  source  only  excites  a  >’-component  of  the  electric  field,  only  a  plot  of  Ey  is 
shown.  Yet  for  distances  as  far  out  as  675  meters,  and  low  vertical  depths,  the  quasi-static 
data  agree  with  the  Sommerfeld  data  reasonably  well. 
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Figure  1:  Electric  field  (y-component)  of  a  2500  A-m:  vertical  magnetic  dipole 
(VMD)  for  a  three  layer,  flat  earth.  The  source  is  in  air  at  a  height  of  15  m  and  the 
observation  point  is  at  a  radial  distance  of  150  m.  The  field  is  a  function  of  vertical 
distance  from  -180  m  in  water  to  45  m  in  air. 
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Figure  2:  Magnetic  field  (z-component)  of  a  2500  A-nT  vertical  magnetic  dipole 
(VMD)  for  a  three  layer,  flat  earth.  The  source  is  in  air  at  a  height  of  15  m  and  the 
observation  point  is  at  a  radial  distance  of  150  m.  The  field  is  a  function  of  vertical 
distance  from  -180  m  in  water  to  45  m  in  air. 
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Figure  3:  Magnetic  field  (x-component)  of  a  2500  A-m2  vertical  magnetic  dipole 
(VMD)  for  a  three  layer,  flat  earth.  The  source  is  in  air  at  a  height  of  15  m  and  the 
observation  point  is  at  a  radial  distance  of  150  m,  The  field  is  a  function  of  vertical 
distance  from  -1 80  m  in  water  to  45  m  in  air. 
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Figure  4:  Electric  field  (y-component)  of  a  45  A-m  horizontal  electric  dipole  (HED) 
for  a  three  layer,  flat  earth.  The  source  is  in  water  at  a  depth  90  m  and  the  observation 
point  is  at  a  radial  distance  of  675  m.  The  field  is  a  function  of  vertical  distance  from 
-180  m  in  water  to  45  m  in  air. 
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Figure  5:  Magnetic  field  (z-component)  of  a  45  A-m  horizontal  electric  dipole  (HED) 
for  a  three  layer,  flat  earth.  The  source  is  in  water  at  a  depth  90  m  and  the  observation 
point  is  at  a  radial  distance  of  675  m.  The  field  is  a  function  of  vertical  distance  from 
- 1 80  m  in  water  to  45  m  in  air. 


Figure  6:  Magnetic  field  (jc-component)  of  a  45  A-m  horizontal  electric  dipole  (HED) 
for  a  three  layer,  flat  earth.  The  source  is  in  water  at  a  depth  90  m  and  the  observation 
point  is  at  a  radial  distance  of  675  m.  The  field  is  a  function  of  vertical  distance  from 
- 1 80  m  in  water  to  45  m  in  air. 
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Data  Comparisons:  Simulation  vs.  Experimentation 


The  plots  on  the  following  pages  show  comparisons  between  data  as  obtained  from 
experimentation  and  simulation.  Actual  experimental  run  numbers  are  shown  in  the  figure 
captions.  Experimental  runs  were  selected  based  on  those  runs  that  had  sufficient 
contiguous  data  streams  unaffected  by  noise  and  random  noise  spikes.  Runs  were 
grouped  to  form  a  single  reference  set  for  comparison.  Due  to  the  way  the  experiments 
were  conducted  and  the  way  the  data  was  collected,  it  was  not  uncommon  to  have  only 
three  data  points  per  set.  Finally,  not  all  of  the  March  2010  experimental  data  was 
processed  during  Phase  11.  Additional  processing  of  that  data  is  also  occurring  in  Phase 
Three;  which  will  be  presented  in  the  Phase  Three  Final  Report.  The  following  table 
correlates  the  figure  numbers  with  the  run  numbers. 


Figures  7,  8  and  9 

Runs  3003,3111  and  3203 

Figures  10,  1 1  and  12 

Runs  2003,  2007  and  2011 

Figures  13,  14  and  15 

Runs  2404,  2408  and  2412 

Figures  16,  17  and  18 

Runs  4304x,  4304y  and  4304z 

Figures  7-9  show  the  electric  field  components  excited  by  a  100  Hz  portable  electric 
source  at  a  depth  of  15.2  meters  at  map  location  4.  The  observation  distance  is  505  meters 
at  map  location  5.  The  data  are  presented  as  a  function  of  sensor  depth.  The  correlation 
between  data  sets  is  quite  good  for  the  Eyand  Ez  components;  the  correlation  is  less  than 
adequate  for  the  Ex  component.  However,  the  modeling  data  for  Ex  are  grouped  together 
with  the  experimental  data  being  the  outlier.  This  suggests  that  the  models  are  consistent 
in  the  way  that  the  experiment  is  interpreted  but  that  interpretation  may  be  wrong.  Further 
study  is  needed  to  assess  and  rectify  this  problem. 

Figures  10-12  show  electric  field  excited  by  a  100  Hz  portable  electric  source  at  a  depth 
of  15.2,  72.2  and  131.7  meters  at  map  location  4.  The  observation  distance  is  964  meters; 
see  map  location  6.  The  data  are  presented  as  a  function  of  source  depth.  In  this  case,  the 
correlation  between  data  sets  is  much  better,  albeit  not  perfect.  The  y-component  of  the 
electric  field  has  two  experimental  data  points  that  lie  near  the  modeling  data;  however, 
the  third  data  point  at  131  meters  is  an  outlier.  The  correlations  for  Ex  and  Ez  are  much 
better.  With  respect  to  Ez,  the  quasi-static  data  is  seen  to  be  off  by  a  factor  of  ten.  Yet, 
since  the  observation  distance  is  964  meters,  such  distances  fall  outside  the  domain  of 
validity  for  the  quasi-static  method. 

Figures  13-15  show  electric  field  excited  by  a  1,000  Hz  portable  electric  source  at  a  depth 
of  15.2,  76.2  and  121.9  meters  at  map  location  4.  The  observation  distance  is  1,000 
meters;  see  map  location  6.  The  data  are  presented  as  a  function  of  source  depth.  The 
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correlation  between  data  sets  is  by  far  the  worst.  However,  the  quasi-static  method 
appears  to  give  the  best  results,  even  though  the  observation  distance  is  far  from  the 
source. 

Finally,  Figure  16-18  show  electric  field  excited  by  the  magnetic  source  on  the  shore  at 
map  location  3.  The  observation  distance  is  675  meters  at  map  location  6.  The  data  are 
presented  as  a  function  of  sensor  depth.  The  HFSS  data  are  clearly  questionable.  This 
poor  data  is  attributed  to  the  way  the  source  is  modeled  in  HFSS;  a  new  source  model  has 
since  been  developed  and  is  currently  being  tested.  The  FDTD  method  seems  to  give 
better  results,  but  the  data  for  Ex  are  particularly  bad.  This  poor  correlation  can  be 
explained  by  noting  that  the  source  is  on  the  shore  and  the  FDTD  method  models  the 
shore  as  if  it  were  a  homogeneous  substance  of  mud,  which  it  is  not.  If  the  constitutive 
composition  of  the  shore  is  not  known  somewhat  precisely  in  the  vicinity  of  the  shore, 
then  there  is  no  expectation  that  the  model  will  predict  the  experimental  outcome. 

To  highlight  this  last  comment,  consider  Figures  19  and  20.  Both  figures  show  the  field 
data,  as  obtained  from  FDTD  simulation,  for  two  different  value  of  shore  conductivity 
(i.e.  rock  vs.  mud);  all  other  parts  of  the  simulation  are  the  same  for  the  two  cases  (i.e. 
source,  geometry,  etc.).  There  is  no  question  from  this  data  that  the  fields  are  highly 
dependent  on  conductivity,  thus  supporting  the  previous  claim  that  when  the  source  (or 
observation  point)  is  near  a  material  boundary  or  interface,  the  constitutive  composition 
of  that  material  must  be  known  to  a  high  degree  of  accuracy  if  good  data  are  to  be 
obtained. 

It  should  be  noted  that  magnetic  field  data  are  not  shown,  even  for  magnetic  source 
excitation.  This  is  due  to  the  very  weak  magnetic  field  signal  that  was  received  relative  to 
the  noise  floor.  Even  when  the  source  and  observation  points  were  close  to  each  other  (i.e. 
100  meters),  the  signal  was  too  weak  to  detect.  The  noise  floor  could  be  reduced  by 
integrating  the  time  domain  data  over  longer  periods  of  time,  but  that  would  require  the 
sensor  array  to  be  stationary  for  long  periods  of  time,  which  it  was  not.  A  Phase  Three 
experiment  will  be  conducted  to  rectify  this  latter  problem  along  with  new 
post-processing  methods  that  will  account  for  sensor  motion. 
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Figure  7:  Electric  field  (^component)  excited  by  a  100  Hz  portable  electric 
source  at  a  depth  of  15.2  meters  at  map  location  4.  The  observation  distance  is 
505  meters  at  map  location  5.  The  data  are  presented  as  a  function  of  sensor  depth 
and  corresponds  to  runs  3003,  3111  and  3203. 


Figure  8:  Electric  field  (x-component)  excited  by  a  100  Hz  portable  electric 
source  at  a  depth  of  15.2  meters  at  map  location  4.  The  observation  distance  is 
505  meters  at  map  location  5.  The  data  are  presented  as  a  function  of  sensor  depth 
and  corresponds  to  runs  3003,  3111  and  3203. 
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Figure  9:  Electric  field  (z-component)  excited  by  a  100  Hz  portable  electric 
source  at  a  depth  of  15.2  meters  at  map  location  4.  The  observation  distance  is 
505  meters  at  map  location  5.  The  data  are  presented  as  a  function  of  sensor  depth 
and  corresponds  to  runs  3003,  3111  and  3203.  The  dotted  line  is  the  experimental 
noise  floor. 
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Figure  10:  Electric  field  (y-component)  excited  by  a  100  Hz  portable  electric 
source  at  a  depth  of  15.2,  72.2  and  131.7  meters  at  map  location  4.  The 
observation  distance  is  964  meters;  see  map  location  6.  The  data  are  presented  as 
a  function  of  source  depth  and  corresponds  to  runs  2003,  2007  and  2011.  The 
dotted  line  is  the  experimental  noise  floor. 
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Figure  11  Electric  field  (x-component)  excited  by  a  100  Hz  portable  electric 
source  at  a  depth  of  15.2,  72.2  and  131.7  meters  at  map  location  4.  The 
observation  distance  is  964  meters;  see  map  location  6.  The  data  are  presented  as 
a  function  of  source  depth  and  corresponds  to  runs  2003,  2007  and  2011.  The 
dotted  line  is  the  experimental  noise  floor. 


Figure  12:  Electric  field  (z-component)  excited  by  a  100  Hz  portable  electric 
source  at  a  depth  of  15.2,  72.2  and  131.7  meters  at  map  location  4.  The 
observation  distance  is  964  meters;  see  map  location  6.  The  data  are  presented  as 
a  function  of  source  depth  and  corresponds  to  runs  2003,  2007  and  2011.  The 
dotted  line  is  the  experimental  noise  floor. 


17 


le-05 


lc-06 


E 

& 


§  le-07 


E 

< 


lc-08 


lc-09 


150 


- 100 

Vertical  distance  (m) 


-50 


Figure  13:  Electric  field  (v-component)  excited  by  a  1,000  Hz  portable  electric 
source  at  a  depth  of  15.2,  76.2  and  121.9  meters  at  map  location  4.  The 
observation  distance  is  1,000  meters;  see  map  location  6.  The  data  are  presented 
as  a  function  of  source  depth  and  corresponds  to  runs  2404,  2408  and  2412.  The 
dotted  line  is  the  experimental  noise  floor. 


Figure  14:  Electric  Field  (^-component)  excited  by  a  1,000  Hz  portable  electric 
source  at  a  depth  of  15.2,  76.2  and  121.9  meters  at  map  location  4.  The 
observation  distance  is  1,000  meters;  see  map  location  6.  The  data  are  presented 
as  a  function  of  source  depth  and  corresponds  to  runs  2404,  2408  and  2412.  The 
dotted  line  is  the  experimental  noise  floor. 
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Figure  15:  Electric  field  (z-component)  excited  by  a  1,000  Hz  portable  electric 
source  at  a  depth  of  15.2,  76.2  and  121.9  meters  at  map  location  4.  The 
observation  distance  is  1,000  meters;  see  map  location  6.  The  data  are  presented 
as  a  function  of  source  depth  and  corresponds  to  runs  2404,  2408  and  2412.  The 
dotted  line  is  the  experimental  noise  floor. 


Figure  16.  Electric  field  (x-component)  excited  by  the  magnetic  source  on  the 
shore  at  map  location  3.  The  observation  distance  is  675  m  at  map  location  6.  The 
data  are  presented  as  a  function  of  sensor  depth  and  corresponds  to  runs  4304x, 
4304y  and  4304z. 
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\e-04 


Figure  17:  Electric  field  (z-component)  excited  by  the  magnetic  source  on  the 
shore  at  map  location  3.  The  observation  distance  is  675  m  at  map  location  6.  The 
data  are  presented  as  a  function  of  sensor  depth  and  corresponds  to  runs  4304x, 
4304y  and  4304z. 


Figure  18:  Electric  field  (y-component)  excited  by  the  magnetic  source  on  the 
shore  at  map  location  3.  The  observation  distance  is  675  m  at  map  location  6.  The 
data  are  presented  as  a  function  of  sensor  depth  and  corresponds  to  runs  4304x, 
4304y  and  4304z. 
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Electric  Field  -  FDTD  -  Case  4  -  Rock  vs  Mud  - 1  KHz 

Rock  ( o^O.OO0 12,  e=4)  Mud  (<M>.01  2,  e=1 ) 


Figure  19:  FDTD  electric  field  data  for  different  conductivities  for  the  shore  and 
lake  bottom  (i.e.  rock  vs.  mud). 
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Magnetic  Field  -  FDTD  -  Case  4  -  Rock  vs  Mud  - 1  KHz 

Rock  (<M).00012,  f=4)  Mud  (o=0.012,  e*l) 
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Vertical  Distance  (m) 


Figure  20:  FDTD  magnetic  field  data  for  different  conductivities  for  the  shore 
and  lake  bottom  (i.e.  rock  vs.  mud). 
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Future  Work 


Based  on  the  previous  discussions  and  the  objectives  of  the  project,  the  following  future 

work  is  envisioned: 

•  Fully  validate  the  new  source  model  in  HFSS. 

•  Implement  a  new  data  post-processing  method  that  accounts  for  sensor  and  source 
rotation. 

•  Conduct  one  more  set  of  experiments  on  the  lake.  Particular  emphasis  will  be  on 
source  and  sensor  locations  near  the  shore  and  on  magnetic  field  sensing  from  the 
magnetic  source  when  source  and  sensor  are  close  to  each  other. 

•  Refine  the  FDTD  and  quasi-static  codes  to  obtain  better  data  as  compared  to  the 
experimental  data. 

•  Improve  FDTD  processing  times  using  parallelization  techniques. 

•  Deliver  Phase  Three  user  manuals  and  documentation  for  each  of  the  developed 
models  and  codes. 

Journals  Articles  to  Date: 

•  Y.  Xia  and  D.M.  Sullivan,  “Underwater  FDTD  simulations  at  extremely  low 
frequencies,”  IEEE  Antennas  and  Wireless  Propagation  Letters ,  vol.  7,  pp. 
661-664,  2008. 

•  Y.  Xia  and  D.  M.  Sullivan,  Z.  Li,  and  R.  Olsen  “Dual  problem  space  FDTD 
simulation  for  underwater  ELF  applications,”  IEEE  Antennas  and  Wireless 
Propagation  Letters,  vol.  8,  pp.  498-501,  2009. 

•  C.  L.  Wagner  and  J.  L.  Young,  “FDTD  numerical  tests  of  the  convolutional-PML 
at  extremely  low  frequencies,”  IEEE  Antennas  and  Wireless  Propagation  Letters, 
vol.  8,  pp.  1398-1401,2009. 

•  D.  M.  Sullivan,  Y.  Xia  and  D.  Butherus,  “A  perfectly  matched  layer  for  lossy 
media  at  extremely  low  frequencies,”  IEEE  Antennas  and  Wireless  Propagation 
Letters ,  vol.  8,  pp.  1080-1083,  2009. 

These  journal  papers  are  appended  to  this  report. 

Papers  in  Preparation: 

•  R.  G.  Olsen  and  Z.  Li,  “A  simple  up-over-down  model  for  low  frequency 
horizontal  electric  dipole  propagation  near  an  interface,”  IEEE  Transactions  on 
Antennas  and  Propagation.  (Appended  to  this  document.) 
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Conference  Publications  to  Date: 


•  D.  M.  Sullivan  and  Y.  Xia,  “Underwater  ELF  simulation  using  the  FDTD 
method,”  IEEE  International  Antennas  and  Propagation  Symposium  and 
USNC/URSI  Radio  Science  Meeting,  San  Diego,  CA,  June  5-1 1,  2008. 

•  Y.  Xia  and  D.  M.  Sullivan,  “Near  to  far  field  transformation  for  underwater  ELF 
simulation,”  IEEE  International  Antennas  and  Propagation  Symposium  and 
USNC/URSI  Radio  Science  Meeting,  San  Diego,  CA,  June  5-11,  2008. 

•  Y.  Xia,  A.  Monsoori,  D.  M.  Sullivan  and  J.  Nadobny,  “High  resolution 
interpolation  for  underwater  FDTD  simulation  at  ELF  frequencies,”  IEEE 
International  Antennas  and  Propagation  Symposium  and  USNC/URSI  Radio 
Science  Meeting,  Charleston,  SC,  June  1-5,  2009. 

•  D.  M.  Sullivan,  Y.  Xia,  “A  perfectly  matched  layer  for  lossy  media  at  extremely 
low  frequencies,”  IEEE  International  Antennas  and  Propagation  Symposium  and 
USNC/URSI  Radio  Science  Meeting,  Charleston,  SC,  June  1-5,  2009. 

•  C.  L.  Wagner  and  J.  L.  Young,  “Characterizing  the  convolutional  perfectly 
matched  layer  at  extremely  low  frequencies,”  2010  IEEE  International 
Symposium  on  Antennas  and  Propagation  and  CNC/USNC/URSI  Radio  Science 
Meeting,  Toronto,  Ontario,  July  2010. 

•  D.  Butherus,  Y.  Xia,  D.  M.  Sullivan,  “Time-domain  near-to-far  field 
transformation  for  underwater  FDTD  simulations  at  ELF  frequencies,”  2010  IEEE 
International  Symposium  on  Antennas  and  Propagation  and  CNC/USNC/URSI 
Radio  Science  Meeting,  Toronto,  Ontario,  July  2010. 

•  Y.  Xia  and  D.  M.  Sullivan,  “Underwater  ELF  Simulation  using  dedicated 
hardware,”  2010  IEEE  International  Symposium  on  Antennas  and  Propagation 
and  CNC/USNC/URSI  Radio  Science  Meeting,  Toronto,  Ontario,  July  2010. 

•  D.  M.  Sullivan,  Y.  Xia,  A.  Mansoori,  “Large  scale  underwater  FDTD  ELF 
simulations  using  Acceleware  and  MPI  parallel  processing,”  URSI  International 
Symposium  on  Electromagnetic  Theory,  Berlin,  Germany,  August  16-19,  2010. 

Master  of  Science  Theses: 

•  Y.  Xia,  “Three-Dimensional  FDTD  Simulation  for  Underwater  ELF  Signals,” 
Master  of  Science  in  Electrical  Engineering,  University  of  Idaho,  August,  2008. 

•  A.  Mansoori,  “Parallel  Three-Dimensional  FDTD  Algorithm  Using  the  MPI 
Library,”  Master  of  Science  in  Electrical  Engineering,  University  of  Idaho, 
August,  2010. 


23 


Appendix  A 

Electromagnetic  Fields  from  an  Electric  or  Magnetic  Dipole  in  a 

Three-layered  Medium 

Robert  G  Olsen  and  Zhi  Li 

School  of  Electrical  Engineering  &  Computer  Science,  Washington  State  University 


Introduction 

The  electric  (E)  and  magnetic  (H)  fields  from  a  dipole  (electric  or  magnetic),  which  is 
placed  in  the  top  half  or  buried  in  the  middle  layer  of  a  three-layer  medium,  can  be 
determined  by  the  Sommerfeld  integral  method.  The  objective  of  this  report  is  to  find  and 
validate  the  solutions  to  the  E  and  H  fields  anywhere  in  the  model  in  terms  of  the 
Sommerfeld  integrals. 

In  the  three-layered  model,  two  half  spaces  occupy  the  top  and  bottom  of  the  medium, 
which  are  denoted  as  #0  and  #2  medium.  The  #0  medium  is  assumed  to  be  free  space  and 
the  #1  is  lossy,  representing  lake  bottom.  And  in  between  there  is  a  layer  conducting 
medium  with  a  uniform  thickness  of  d.  The  middle  layer  is  denoted  as  #1  medium  and  it 
represents  lake  water.  In  this  project,  the  dipole  source  is  allowed  to  be  placed  in  either 
#0  or  #1  medium.  Therefore,  according  to  the  type  (electric  or  magnetic),  orientation 
(vertical  or  horizontal)  and  position  (in  #0  or  #1  medium),  there  are  eight  different  cases 
of  dipole  source  to  be  studied  in  this  project.  The  eight  cases  and  their  assigned  identifiers 
are  shown  in  Table  1 . 

Table  1  Eight  cases  of  dipole  source 


Case 

number 

Descriptions 

Identifier 

i 

Horizontal  electric  dipole  (HED)  in  #0  medium; 

HED0 

2 

Horizontal  electric  dipole  (HED)  in  #1  medium; 

HED1 

3 

Horizontal  magnetic  dipole  (HMD)  in  #0  medium; 

HMD0 

4 

Horizontal  magnetic  dipole  (HMD)  in  #1  medium; 

HMD1 

5 

Vertical  electric  dipole  (VED)  in  #0  medium; 

VED0 

6 

Vertical  electric  dipole  (VED)  in  #1  medium; 

VED1 

7 

Vertical  magnetic  dipole  (VMD)  in  #0  medium; 

VMD0 

8 

Vertical  magnetic  dipole  (VMD)  in  #1  medium. 

VMD1 
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Without  losing  generality,  the  case  of  HMD  in  #1  medium,  HMD1  is  chosen  as  example 
in  the  report  to  show  the  detailed  process  how  to  find  the  solutions,  by  the  Sommerfeld 
integral  method,  to  E  and  H  fields  at  the  field  point  anywhere  in  the  model.  The  results  of 
the  other  seven  cases  will  be  listed  at  the  end  of  the  report. 

Geometry  of  Case  HMD1 

For  the  case  HMD1,  a  horizontal  magnetic  dipole  (HMD)  is  placed  in  the  middle  layer 
(#1)  of  the  three-layered  medium.  The  ly*  oriented  HMD,  with  a  dipole  moment  of  IdA 
(A-nT),  is  on  the  ‘z’  axis  and  buried  in  medium  #1  and  ‘ h ’  meters  below  the  interface 
between  medium  #0  and  #1  (- d  <  z  <  0).  The  cylindrical  coordinate  system  (p,  <p,  z)  is 
used  in  this  paper,  where  x  =  pcosy  and  y  =  ps\n(p.  Thus,  the  observation  (or  field)  point  is 
assumed  to  be  at  (p ,  (p,  z).  Fig.  1  illustrates  the  model. 


Free  Space 
(#0) 

,z 

z  =  0 

«o.  Ob.  M, 

y 

observation 

« 

HMD=*Z  =  -h  *  P°'nt 


Conducting 
Medium  (#1) 

. cl  >  O) ,  //0 

777777777 

^77^77777777777 ~ 

Bottom 

C2>  02,  //0 

(#2) 

Fig.  1  Illustration  of  the  model 

As  noted  in  the  figure,  e,  and  a,  are  the  permittivity  and  conductivity  of  the  ilh  half  space 
(/  =  0  and  1  for  free  space  and  conductor,  respectively),  e,  =  cne o,  where  is  the  relative 
permittivity  and  eo  is  the  permittivity  of  free  space.  It  is  assumed  that  all  materials  have 
the  permeability  of  free  space  po-  Before  the  derivations,  some  of  the  constants  and 
variables  need  to  be  defined. 

=  -co2pe]  =  -«>(£•,  -j—) 

CO 

u,=Ju2+r?) 

where  e\-=e,- joja  is  the  complex  permittivity  of  the  medium  #/’ s  (/  =  0,  1  or  2),  k,  is 

the  wave  number  where  Re(y,)^0  and  Re(w,)  SO  defines  the  proper  Reimann  sheet  of  the 
complex  plane. 
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Sommerfeld  Integral  Equations  to  E  and  H  Fields 

According  to  the  theory  of  electromagnetics,  the  fields  due  to  the  dipole  radiation  should 
satisfy  the  boundary  condition  at  both  the  interfaces.  For  the  horizontal  dipole  case,  to 
fulfill  the  boundary  conditions  two  non-zero  components  of  vector  potentials  are  required. 
The  pair  can  be  a  variety  of  combinations  of  the  vector  potentials.  In  this  project,  for  the 
convenience  of  derivation,  we  choose  Fy  and  Fz  to  be  the  only  non-zero  vector  potentials. 
The  vector  potentials  can  be  written  in  terms  of  Sommerfeld  integrals  as: 

F;  =  kx^MX)e-u-XJ0(Xp)dX  (z  >  0)  (1) 

FhK^  +  K  f[/2(A)<f +  MW''- ] U0(Ap)dA  (~d<z<  0)  (2) 


Fy  =  *,  £  (z  <  -d)  o) 

enR  [  £  u;'e-*'u+h)M0(Xp)dA  (z  +  h)  >  0 

where  - =  <  is  the  source  term  and  represents  the 

R  J ~ux'eu'(z*h)XJ0{Xp)dX  (z  +  h)<0 


field  directly  from  the  dipole  source,  R  =  (p2  +  z2)1'2  is  the  distance  from  the  dipole  to  the 


observation  point,  and  The  superscripts,  0,  1  or  2,  of  vector  potential 

An 


indicates  w  hich  layer  of  medium  the  vector  potential  is  related  to. 

Fz  =  f  (z  -  0) 


(4) 


F!  [[g.We-*'*  +g,(\)eu'']XJQ(Xp)dA  (~d<z<  0)  (5) 

F2  =  ^  f (z  <  -d)  (6) 

Since  the  dipole  is  oriented  in  'y*  direction,  only  the  F'  component  contains  source 

term,  as  shown  in  (2).  In  these  equations,/;  ~  f4  and  gi  ~  g4  are  just  some  arbitrary 
coefficient  functions  of  the  integral  variable  X  to  be  determined  by  the  boundary 
conditions.  To  find  the  solutions  to  the  fields,  the  first  step  is  to  determine  these 
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coefficient  functions. 


Coefficient  Functions 

Known  the  vector  potentials,  the  electric  and  magnetic  fields  can  be  calculated  by 
Maxwell’s  equations  (7)  and  (8). 


Ef  =--VxF 
£ 


(7) 


HF  =  -jcoF - — V(V  •  F) 

cofue 


(8) 


Equation  (9)  to  (14)  give  the  expressions  of  each  field  components  in  terms  of  Fy  and  Fz. 


e  dy  dz 


r,  _  1  dF, 
y  e  dx 


1  dF 
F  =  y 

e  dx 


//.  — 


j  d  (dFL  +  dF^) 


(Ofis  dx  dy  dz 


COflS 


dy  dy  dz 


H.=  — 


CO/.IS 


'  d  dF  dF  2„ 
— ( — -  +  — -)-r  F. 
dz  dy  dz 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 


The  boundary  conditions  to  be  satisfied  are  that  all  the  tangential  fields  are  continuous  at 
both  the  upper  and  lower  interfaces.  At  the  upper  interface  (z  =  0),  they  can  be  written  as: 


(15) 
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(16) 


i  bf;  =  i  8f; 

e'0  dz  e[  dz 


1 

f  5F>°  ,  3F2° ) 

1 

f  dFy  dF!) 
— -  +  — 

*0 

dz  j 

£'1 

{  dy  dz  J 

y_ o  r^o y 1  77  * 

£•'  '  £'  y 

*0  cl 

And  at  the  lower  interface  (z  -  -ii),  they  are: 


1  dF[  _  1  dF? 
e'2  dz  e[  dz 


dF :  dF.2 


dy  dz 


dF'  dF' 


dy  dz 


(17) 

(18) 

(19) 

(20) 

(21) 


rLF2  =  rLF' 

s'  y  s'  y 


(22) 


Plugging  all  the  vector  potentials  into  (15)  to  (22)  and  simplifying  the  equations,  it  is 
obtained  that 


£“0  £x 


(23) 


^  f\  W  =  -7  [e~",h  +  w./2  (A)  ~  «,/,  (A)]  (24) 

£*0  £\ 


-tI/w-w-oH-t 

£*n 


— e 


+  fi  (A)  +  f}  (A)  -  z/,£2  (A)  +  uxg3  (A) 


(25) 


yj(A)  =  —e~u'h  + /2(A)  4-/3  (A) 


(26) 
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\e-u'dg^)  =  ±[eu'dg2{X)  +  e-u'dg3{*.)] 


=  -L[e“'<W) 


(27) 


(28) 


—Utd 

e  1 


— [/4(A)+w2g4(A)]  =  — 


eUidf,  (A)  =  +  e*f2  (A)  +  e-'/a  (A) 

u, 


eu'  '  +  eU|  /2  ( A)  +  e_U|  /3  ( A)  -  M,eU|  g2  ( A)  +  g3  (A) 

(29) 

(30) 


Let  (24)  -  (26)  x  — 2-  to  cancel  the  term  of //,  we  have 


t  t 

f2  + 

t  \ 

$  f 

fl=- 

(  \ 

Ml_Ml 

t  f 

V  s\  £o  J 

\  £\  £o) 

U\ 

V  *0  £\  ) 

-uxh 


Similarly,  let (28)- (30) x-j-  to  get  rid  of /)  such  that 

S 


u,  u-, 


£\  £2  y 


e-"7,+R-T  e-*7,-LR-* 


f,  f2  y 


Mi  v  ^2  y 


u,(W) 


(31) 


(32) 


Solve  equation  (31)  and  (32)  together  for  /  and  /j 

fi  =  -^\i.^-^oui){e2u\-e'iui )e’“'(W)  +«w0  +^«l)(e1'M2-^M1)eu,(W)]  (33) 


/3  =^[(£0M1  -ff,'«o)(ei'w2 -^miK'<W)+(S|mo -£omi)(^'u2  +^M,)e'“l<W)]  (34) 
where  Z)  =  (e'l^ )(£■'«, -£l,2/2)e_U|‘/ -(f'w,  +  £1'w0)(£,'w2  +  ) eU|</ .  Then  use  (26) 


and  (30)  to  find//  and /) 


2s'eu'J 


f*  =  2C—[{s\-s'au,)e-u'h  -(e'0u ,  +  s,iu0)eu'h'\ 


(36) 


Next  step  is  to  solve  equations  (27)  to  (30)  for  gi  ~  g4.  Before  that,  it’s  convenient  to 
simplify  (28)  with  (24) 


i  ,  v 


*0 

T£i  = 


r  w.  u, 
f\  — 7^2  +,£3 


£\  £\ 


(37) 


and  to  simplify  (30)  with  (26) 


<2  -11,6/ 

Te  "  £4  = 


Ve>  ^2  7 


r  u.c/  .  “l  -u.c/ 

•/. — re  £2  +  1  £3 


(38) 


It  is  the  similar  routine  to  getg^  and  gs,  then  gi  and  g4.  Let  (5 )xw0  +(15) : 

\ 

fx 


(s’  ^ 

0<O-Ml)£2+(MO+Wl)£3  =  — 


y 


(39) 


and  (7)xh2+(16): 


(«2  +  «,)*"''  '£2  +  («2  ‘£3  =1  l-^ 


/4 


'2  / 


Solve  (39)  and  (40),  we  have 


£2  = 


_  g;  (fj -  go)(^2  -  »i  )g"V  f\  +  go(gl'  -  g2  XM0  +  »i  K"1  /4 


£0£2  •  A 


_  ^0  ~gl')(»l  +  -/j  +go(g2  -gl')(»0  Z3K1  li 

£0£2  DX 


(40) 


(41) 


(42) 


where  Dt  =  (w0 -mi)(m2  -ut)e  -(u0  +  ul)(u{+u2)eu,‘ .  Then  solve  (23)  and  (27)  to 
obtain  gi  and  g4. 

_  g'(g;-g;)[(»2  -(2/,  +u2)e’*l]fl +2g>,(g;-g;)g-^/4 


£.  =• 


£,A  A 


(43) 
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(44) 


8  4  = 


_  2g>|  (gp  -  +  gp  (f [  -  g;  )  [(«0  +  »I  )gl'l</  -  (»0  -  MI  )g"''l</  ]  h 


£0£{ 


D, 


Integral  Equations  of  E  and  H  Fields 

Since  all  the  coefficient  functions  have  been  determined  in  part  2.1,  the  E  and  H  fields 
anywhere  in  the  model  can  be  determined.  After  plugging  the  vector  potentials  in  (1)  to  (6) 
into  the  Maxwell's  equations  (9)  to  (14)  and  some  algebra  manipulations,  the  components 
of  E  and  H  fields  at  the  observation  point  can  be  obtained.  The  results  are  shown  below. 
In  medium  #0  and  #2,  there  is  no  dipole  source.  Therefore,  the  field  component  contains 
only  the  transmitted  field  from  the  interface.  The  transmitted  field  is  denoted  by  a  7’  in 
subscript  of  each  field  component. 

The  field  component  in  medium  #1  is  formed  by  two  parts:  one  is  the  incident  field 
directly  from  the  dipole  source  and  the  other  one  is  the  reflected  field  due  to  the  two 
interfaces.  They  are  denoted  by  7’  and  V’  in  the  subscript  of  the  corresponding  field 
component. 
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In  medium  #0  the  results  are: 
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In  medium  #1  the  integral  equations  for  the  E  and  H  fields  are: 
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Numerical  Validation  of  the  Results 


To  validate  the  results  in  (45)  to  (62),  the  electric  and  magnetic  are  numerically  calculated 
and  checked  with  some  known  results,  such  as  the  quasi-static  fields  in  infinite  medium. 
The  values  of  the  parameters  used  in  the  following  calculation  are  listed  in  Table  2. 


Table  2  List  of  parameters 


Medium 

#0 

#1 

#2 

Relative  permittivity,  sr 

1 

81 

3 

Conductivity,  a  (S/m) 

0 

4 

0.01 

Permeability,  n  (H/m) 

4tcx10-/ 

4:1*1 0‘v 

471*1  O'7 

d(m) 

300 

h(m) 

can  vary  from  0  to  300 

Dipole  moment  IdA  ( A  •  m2 ) 

i 

Frequency/ (Hz) 

10  to  3000 

Numerical  Integration 

In  this  project  we  choose  the  composite  Simpson's  rule  to  do  the  numerical  integration. 
For  a  given  integrand /(6r),  the  Simpson’s  rule  is  used  to  obtain  the  integration  of f(x)  over 
interval  [a,  b\.  It  is  given  as 

f  mdx = |  [  m +4 ■/(*,)+ /(*)]  -^  /(4)  (4)  (63) 

where  xi  is  the  middle  point  of  [a,  6],  a<£  <b .  The  Simpson’s  rule  is  usually  inaccurate 
if  used  over  large  integration  intervals.  To  avoid  the  problem,  a  piecewise  approach,  the 
composite  Simpson’s  rule,  is  often  applied.  (Fig.2) 
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First,  the  integration  interval  [a,  £]  is  divided  into  n  equal-spaced  subintervals,  where  n 
must  be  an  even  number.  Then  apply  the  Simpson’s  rule  on  each  subinterval  and  combine 
all  the  integrations  over  every  subinterval  to  get  the  final  result.  In  formula,  the  composite 
Simpson’s  rule  is  described  as 


in!  2H 


nil 


/W+2  I  /<^)+4l/(^H)+/(JI,) 

>-l  J- 1 


]P  "/2 

<«) 


y-i 


where  x2j_2  <,  <,  x2J ,  for  each  j  =  1,  2,  n/2.  When  the  numerical  integration  is 

carried  out,  the  error  term  is  usually  truncated. 


(«/ 2H 


nil 


/(*„)+ 2  E  /(V+4E /(*=/-.) +/w 

j=i  j.i 


(65) 


Theoretically,  as  shown  in  (45)  to  (62),  the  exact  fields  will  be  given  by  the  integration 
from  zero  to  infinity.  But  it  is  not  possible  to  do  this  in  a  numerical  manner.  The  computer 
program  can  only  deal  with  integration  over  finite  intervals.  To  make  the  fields 
calculation  possible,  some  approximation  should  be  made.  An  integral  can  be  separated 
into  two  parts 

£  f(x)dx  =  £  f(x)dx  +  £  f(x)dx  (66) 


If  we  can  find  a  bound  number  'b ’  such  that  the  second  integral  on  the  right  hand  side  is 
small  enough  compared  to  the  first  integral,  the  total  integral  can  be  approximated  by  the 
first  term.  Fortunately,  we  do  can  find  such  kind  of ‘A’  for  the  field  calculation  because  all 
the  integrands  in  the  field  calculation  equations  have  attenuation  characteristics. 

frX  max  fee 

Fl(X)dX  =  f  FI(X)dX+\  FI(X)dX 

J)  •'/l  max  /  rn\ 

,4-nax  (67) 

*  \  FI{X)dX 

where  FI(X)  represents  integrand  for  field  integration,  Xmax  is  the  upper  limit  of  the 
integration  interval  to  be  used  in  numerical  calculation. 


In  practice,  we  use  50  as  the  value  of2max  for  all  the  terms  derived  from  the  source  term, 
er,*  f  f  u;'e-u'l”h)AJ0(Ap)dA  (z  +  h)  >  0 

R  £  u;'eu'iz+h)AJ0(Xp)dX  (z  +  h)  <  0 

This  value  of  Xmax  will  give  us  enough  accuracy  for  calculation.  However,  for  all  the 
scattered-field  terms,  which  contains  the  coefficient  functions  fi~f4  and  gi  ~  g4 ,  the  same 
An, ax  doesn’t  work.  If  the  2max  value  is  too  large,  the  calculation  of  the  coefficient  functions 
will  exceed  the  operation  limit  of  the  computer  (like  103‘4).  By  testing,  the  Xmax  is  set  at 
1.2.  Although  it  is  much  smaller  than  that  for  the  source  terms,  the  final  results  of 
calculation  are  still  acceptable. 
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The  integral  step  h  is  another  important  factor  to  the  numerical  integration.  Too  large 
steps  not  only  bring  big  error  but  also  cause  bad  behavior  of  the  calculation.  Small  steps, 
however,  slow  down  the  speed  of  integration.  In  this  project  we  choose  the  integral  step  h 
around  0.01  for  both  source  term  and  scattered  term  integration.  And  it  gives  a  good 
compromise  between  computing  stability  and  speed. 

If  the  HMD  source  is  located  in  the  middle  of  the  water  layer  and  we  consider  the  area 
not  too  far  away  from  the  dipole,  the  scattered  fields  in  this  area  will  be  so  small  that  can 
be  ignored.  The  total  field  then  will  behave  exactly  like  that  induced  by  a  HMD  radiating 
in  an  infinite  uniform  media  of  water.  In  order  to  verify  the  numerical  results,  those 
results  are  compared  to  the  fields  radiated  by  a  HMD  in  the  infinite  water  media. 

Fig.  3  and  Fig.  4  show  an  example  of  the  comparison  between  the  two  set  of  results.  The 
HMD  is  at  h  =  150m,  which  is  in  the  middle  of  the  water  layer.  The  ‘0’  angle  of  the 
evaluating  points  is  tt/4  and  z  =  -149m.  For  a  HMD  radiating  in  the  infinite  uniform 
conducting  medium,  Ey  component  is  always  zero.  So  there  is  no  Ey  presented  in  the 
comparison. 

In  those  figures  the  solid  line  curve  stands  for  the  3-layer  results  and  the  star-line  curve 
for  the  uniform  media  results.  It  is  clear  that  the  two  sets  of  the  results  match  each  other 
very  well. 
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h  ■  150>  ■  025x 


h*l«**Q35t 


7 


(a)  E„ 


(b)  Ez 


Fig.  3  Comparison  of  E  fields  between  3-layer  results  and  uniform-media  results 


MUO.fOtti  fi  ■  ISO.#  ■  0  25* 


(a)  H*  (b)  Hy 


(c)  Hz 


Fig.  4  Comparison  of  H  fields  between  3-layer  results  and  uniform-medium  results 
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Checking  Boundary  Conditions 


Checking  boundary  conditions  at  each  interface  provides  another  approach  to  validate  the 
the  solutions  to  the  fields,  equation  (45)  ~  (62).  At  the  upper  interface  (z  =  0),  the 
tangential-field  components  are 

K  r  /•  -  1  AJ0(Ap)sin2(^)  +— J,(/i/?)cos(2^) 

P 


K  =  — r  £  uJ^^p)dX+-L  £  g, 


A2dA 


4iw*> 


/l./0(A/>)sin2(^)+-Jl(A/>)cos(2^) 
P 


A2dA 


£. 


AJ0{Ap) - J^(Ap) 

P 


A2  sin(^)cos (<f>)dA 


£fW*>) 


AJ0{Ap) - J^Ap) 

P 


A 7  s\n(<j>)cos(<j))dA 


AJ0(Ap) - J,(Ap) 

P 


A2  sin(^)cos(^)t//l 


(u^e~u'h  +  fi  +  A  ~ «,g2  +  ".#3 ) 


AJ0(Ap) - J,(Ap) 

P 


A2  sin(^)cos(^)t/A 


AJ0(Ap)  sin2  (</>)  h —  J,  ( Ap)  cos(2  <j>) 


P 


A  dA 


IdA  s! 


An  e 


^^yIWMpWa 


=  (u;'e~y'h  +  f2  +  f3  -  utg2  +  ulgi ) 


AJ0(Ap )  sin2  ( <f>)  + — J,  (Ap)  cos(2^) 
P 


A2dA 


IdA 

An 


[(u;'e~u,h  +f2+f})-rfAJ0(Ap)dA 


With  the  boundary  conditions  (23)  ~  (26),  it  is  not  difficult  to  prove  that 


E°x  =  E[  H°x  =  H'x 

and 

El  =  E[  H°  =  Hl 
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So  the  boundary  conditions  are  theoretically  proved  to  be  satisfied  at  the  upper  interface. 
Similarly,  at  the  lower  interface  (z  =  -d)  the  tangential-field  components  are 

E;  =  4  f  u2e~^fAXJ0(Xp)dX+^  JT  e^g4  ■ 


XJ0(Xp)sin"(4>)  +—Ji(Xp)  cos(2^) 
P 


X2dX 


e,  A  f  [e'M  +u]e~*f,]U„aP)dA 


A  f  («'"*!  +*^*1) 


XJ0  (Xp)  sin 2  (<f>)  +  —  Ji  ( Xp )  cos(2  <j>) 


X2dX 


*2 


XJ0(Xp) - Jx(Xp) 

p 


X?  sin(^)cos (<j>)dX. 


XJQ  (Xp)  ,/|  ( Xp) 

p 


X2  sin(^)cos(^)c/A 


XJ^P) — AUp) 
p 


X2  sin(^)cos(^)c/A 


Hi  =  ^  f  +  e'V;  +«"'7,  -  ) 


ajMp) — 

p 


X2  sin(^)cos (<j>)dX 


XJ q  (/fy?)  sin 2  (^)  +  —  J,  (A/?)  cos(2^) 
P 


X2dX 


IdA  s'. 


4  n  s 


4f  rie-'dfAXJ0(Ap)dX 


H\  = ~  £  (»r'^(W)  +  e*f2  +  <f"'73  - ufi«g2  +  ) 


4;r 


IdA 

An 


XJ0  (X.p)  sin2  (<f>)  + — J,  (Xp)  cos(2^) 


X2dX 


£  («r'^,(W)  +  ^72  +  e^df3)-y2XJ0(Xp)dX 
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By  considering  the  equations  (27)  ~  (30),  it  also  can  be  proved  that 


Hi  =  H\ 


"1  =  H'r 


Then  the  boundary  conditions  are  theoretically  proved  to  be  satisfied  at  the  lower 
interface.  Fig.  5  and  Fig.  6  show  the  matching  of  boundary  conditions  at  the  lower 
interface.  The  HMD  source  is  put  at  h  =  290m.  For  the  field  points,  0=  7r/4.  (Figures  are 
on  next  page.)  As  shown  in  the  figures,  the  fields  at  the  both  sides  match  each  other 
very  well. 
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H»29C.»»025« 


»>SCft2i 


[id ' 


(a)  Ex 


(b)  Ey 


Fig.  5  Checking  tangential  E  fields  along  the  water-bottom  interface 


(a)  Hx 


(b)  Hy 


Fig.  6  Checking  tangential  H  fields  along  the  water-bottom  interface 
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Appendix: 


Sommerfeld  integral  equations  of  E  and  H  fields  for  other  dipole  cases 

There  are  eight  cases  of  different  dipole  source  conditions  interested  in  this  project.  In  the 
previous  sections  of  the  report,  the  process  of  finding  the  Sommerfeld  integral  solutions 
to  E  and  H  fields  for  the  case  that  HMD  in  #1  medium,  HMD1,  has  been  introduced.  By 
using  the  similar  method,  the  solutions  for  the  other  seven  cases  can  also  be  obtained.  In 
this  appendix,  the  results  of  the  cases  (in  case  identifier  as  shown  in  Table  1):  HED  , 
HED1,  HMD0,  VED°,  VED1,  VMD°,  and  VMD1  are  provided. 


45 


HED  -  Dipole  in  #0  medium  (free  space) 
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HED  -  Dipole  in  #0  medium  (free  space) 
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HED°  --  Dipole  in  #0  medium  (free  space) 
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Appendix  B 

Maximum  Detectable  Range  for  Electromagnetic  Fields  from 
Dipole  Sources  Near  an  Air- Water  Interface 

Robert  G.  Olsen  and  Zhi  Li 

School  of  Electrical  Engineering  &  Computer  Science,  Washington  State  University 

Introduction 

In  Phase  II  of  this  Project,  a  Matlab  program  based  on  the  Sommerfeld  integral 
formulation  to  calculate  the  electromagnetic  fields  from  electric  or  magnetic  dipole 
sources  in  a  three-layered  medium  was  developed.  The  program  allows  four 
types/orientations  of  dipole  sources  which  are:  vertical  electric  dipole  (VED), 
horizontal  electric  dipole  (HED),  vertical  (VMD),  and  horizontal  (HMD)  magnetic 
dipole.  The  three  layers  are  numbered,  from  the  top  to  the  bottom,  as  layer  0,  1,  and  2, 
respectively.  The  top  (0)  and  the  bottom  (2)  layers  extent  to  +/-  oo  respectively.  Layer 
0  is  assumed  to  be  free  space,  while  layers  1  and  2  are  conducting  media  and  the 
dipole  source  can  be  placed  in  either  layer  0  or  1.  Fig.  1  illustrates  the  case  for  an 
HED  in  layer  1. 


Free  Space 
(#0) 

,z 

z  =  0 

£o-  °b>  Vo 

y 

observation 
•  •  . 

HED4z=-h  *  P°'nt 


Conducting 

Medium  (#1) 

Si.  ®i.  Vo 

777777777 

Z  v.d  ^77777777777 

Bottom 

c2,  o2,  m, 

m 

Fig  1.  Model  of  a  HED  placed  in  #1  medium. 


Detectable  Range 

Using  this  program,  the  electric  (E)  and  the  magnetic  (H)  fields  anywhere  in  space 
can  be  calculated.  This  provides  the  basis  for  determining  the  detectable  range  from 
the  source  if  the  maximum  dipole  moment  and  the  minimum  detectable  signal  for  the 
measuring  equipment  are  given.  In  these  simulations,  it  was  assumed  that  the 
maximum  dipole  moments  for  electric  and  magnetic  dipoles  are  50  A-m  and  2500 
A-nr  respectively  and  that  the  minimum  detectable  electric  and  magnetic  fields  are 
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lpV/m  and  40pA/m  respectively.  It  was  observed  from  the  simulations  that  the 
horizontal  electric  field  component  perpendicular  to  the  HED  direction  can  nearly 
always  be  detected  at  a  distance  much  larger  than  that  for  any  other  field  component 
from  any  other  dipole  source1. 


Figs.  2-10  show  the  variation  of  the  maximum  detectable  range  of  Ex  with  dipole 
frequency.  Figs.  2  -  4  are  the  results  for  the  cases  that  the  dipole  is  2  meters  below  the 
upper  interface  and  the  field  point  is  5,  20  and  50  meters  below  the  upper  interface, 
respectively.  The  dipole  moment  for  these  three  simulations  is  Idl  =  50  A-m. 


h  =  2;  z  =  -5 


Fig.  2.  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  (HED  is  2  meters  below  the  upper 
interface,  h  =  2,  field  point  is  5  meters  below  the  interface,  z  -  -5). 


1  For  a  conductivity  of  4  8  S/m,  the  HED  magnetic  field  has  a  higher  detectable  range  for  frequencies 
less  than  500  Hz.  At  smaller  conductivities,  the  frequency  at  which  this  occurs  is  smaller  than  this. 
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h  =  2;  z  =  -20 


Fig.  3.  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  (HED  is  2  meters  below  the  upper 
interface,  h  =  2,  field  point  is  20  meters  below  the  interface,  j  =  -20). 


h  =  2;  z  =  -50 


Fig  4  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  (HED  is  2  meters  below  the  upper 
interface,  h  =  2,  field  point  is  50  meters  below  the  interface,  z  ~  -50). 
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In  Figs.  2  and  3,  the  curves  at  the  top  and  bottom  represent  the  simulation  results  for 
lake  water  and  sea  water.  The  conductivity  for  the  curve  in  the  middle  is  an  arbitrarily 
chosen  value  between  that  for  lake  water  and  sea  water.  The  detectable  range  for  sea 
water  when  the  field  point  is  at  z  =  -50m  is  so  short  that  it  is  not  shown  in  Fig.  4. 

The  simulations  for  Figs.  5  to  7  are  respectively  similar  to  that  for  Figs.  2-4.  The 
only  difference  is  that  the  dipole  moment  and  the  minimum  detectable  field  are 
assumed  to  be  increased  and  reduced  by  three  times,  respectively. 

Idrw  =  3  x  Idl 
Enew  =E  /3 

mtn -detectable  min -detectable 

The  detectable  ranges  in  these  cases  are  much  larger  than  their  counterparts  in  Figs.  2 
-  4  due  to  the  increase  of  both  the  strength  of  source  signal  (i.e.,  the  dipole  moment) 
and  the  ability  of  detection  (i.e.,  minimum  detectable  field).  In  Figs.  8  -  10,  the  dipole 
moment  and  ability  of  detection  are  increased  by  five  times: 

Wr=5xM 

Enew  =E  /  5 

min -detectable  mtn  -detectable 

The  increases  of  the  detectable  range  of  Ex  become  more  obvious. 


h  =  2;  Z  =  -5 


Fig.  5.  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  when  the  dipole  moment  and  the 
detectability  are  both  increased  by  three  times  (Mnew  =  3 *Idly  Eminnew  =  EmJ3;  h  -  2,  and  z  “  -5). 
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h  =  2;  z  =  -20 


Fig.  6.  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  when  the  dipole  moment  and  the 
detectability  are  both  increased  by  three  times  (Idlne w  =  3 *Idl,  Emin>new  =  EmJ 3;  h  =  2,  and  z  =  -20). 


h  =  2;  z  =  -50 


Fig.  7.  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  when  the  dipole  moment  and  the 
detectability  are  both  increased  by  three  times  (Idl^  =  3 *Idl,  Emmrww  =  Emin/3 ;  h  =  2,  and  r  =  -50). 
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h  =  2;  z  =  -5 


Fig.  8.  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  when  the  dipole  moment  and  the 
detectability  are  both  increased  by  five  times  (Idlnc w  =  5 *Idl,  EmmnfW  =  EmJ 5;  h  =  2,  and  r  =  -5). 


h  =  2;  z  =  -20 


Fig  9.  Detectable  range  of  as  a  function  of  dipole  frequency  when  the  dipole  moment  and  the 
detectability  are  both  increased  by  five  times  ( IdIMW  -  5 *Idly  Eminn€W  -  EmJ 5;  h  -  2,  and  z-  -20). 
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Fig.  10  Detectable  range  of  Ex  as  a  function  of  dipole  frequency  when  the  dipole  moment  and  the 
detectability  are  both  increased  by  Five  times  (Idlne w  =  5*Idl,  Eminnew  =  EmJ 5;  h  =  2,  and  z  =  -50). 


llp-over-and-down  Model 

Another  part  of  work  accomplished  in  Phase  II  is  the  study  of  an  up-over-and-down 
model  for  HED  propagation  near  an  interface  [1].  In  this  study,  the  HED  is  assumed 
to  be  buried  in  the  low  er  (lossy)  half  medium  of  a  two-half-space  model.  Note  that  the 
bottom  layer  of  Fig.  1  can  often  be  neglected  if  source  and  field  point  are  much  closer 
to  the  top  interface  than  the  lower  one.  The  top  half  is  free  space.  The  model  is  shown 
in  Fig.  1 1. 


Free  space  (#0) 

\Z 

z  =  0 

Co.  Mi 

y 

observation 

HED  - 

»  z  =  -h 

point 

Conducting  medium  (#1) 

i 

£i>  °i>  fh 

Fig.  II. 

Geometry  of  the  model 

When  the  HED  and  the  field  point,  (pt  cp,  z),  are  both  close  to  the  interface,  the 
Sommerfeld  integrals  for  the  electric  and  magnetic  fields  in  the  lower  half  space  can 
be  simplified  and  a  set  of  simple  approximations  for  the  fields  obtained.  For  example, 
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if  the  depths  of  the  HED  and  field  point  are  much  smaller  than  the  horizontal  distance 
between  them,  the  ^-component  of  electric  field  in  the  conducting  medium  can  be 
approximated  by 


E\*2jAxk\-e 


A  (--*). 


-A  p 


3  +  - 


Ukp)  i  0k\p) . 


sin  <f>  cos  (j> 


(1) 


where  a.  =  ^  ■  Equation  (1)  was  used  to  do  the  similar  simulations  for  determining 
1  4  ncoe\ 

the  detectable  range.  Fig.  12  and  Fig.  13  show  the  comparisons  between  the  results 
obtained  by  (1)  and  that  found  by  Sommerfeld  integral  method,  shown  in  Fig.2  and 
Fig.  5,  respectively.  It  is  obvious  that  the  approximation  of  the  field  in  (1)  gives  very 
good  result  over  the  most  portion  of  the  frequency  range  interested.  Since  the 
approximation  (1)  has  no  integral  in  it,  the  calculation  time  can  be  significantly 
reduced  by  using  it.  Therefore,  the  approximation  based  on  the  up-over-and-down 
model  provides  us  a  fast  but  relatively  accurate  approach  to  determine  the 
electromagnetic  field  in  the  conducting  medium  when  the  third  layer  can  be 
neglected. 
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Fig  12.  Comparison  between  Sommerfeld  method  and  up-over-and-down  approximation 
for  detectable  range  of  Ex.  (HED  is  2  meters  below  the  upper  interface,  h  =  2,  field  point  is  5 
meters  below  the  interface,  z  =  -5). 
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Fig.  13.  Comparison  between  Sommerfeld  method  and  up-over-and-down  approximation 
for  detectable  range  of  Ex.  The  dipole  moment  and  the  detectability  are  both  increased  by 
three  times  (/t^/ncw  3 */c//,  Emmnew  Em,J 3,  h  2,  and  z  -5). 
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Appendix  C 

FDTD  Parallelization  Methods 

Dennis  Sullivan 

Department  of  Electrical  and  Computer  Engineering,  University  of  Idaho 

Introduction 

In  the  past  year,  the  FDTD  simulation  group  has  accomplished  the  following: 

1.  The  FDTD  codes  were  reformulated  using  the  Message  Passing 
Interface  (MPI)  software.  The  purpose  was  to  achieve  more 
efficient  parallelization  of  the  code.  It  was  also  to  make  the 
code  more  flexible  by  allowing  it  to  be  distributed  over  several 
computers  [1,2]. 

2.  The  FDTD  codes  were  implemented  using  dedicated  hardware 
from  the  Acceleware  corporation.  This  system  uses  graphics 
cards  to  do  the  bulk  of  the  FDTD  calculation.  This  is  controlled 
through  the  software  development  kit  (SDK)  from  Acceleware  [2, 

3], 

3.  The  near-to-far  field  formulation  was  implemented  in  the  time 
domain  instead  of  the  frequency  domain.  This  method  allows 
greater  flexibility  and  provides  more  information.  Wavelet 
theory  was  used  for  data  compression  to  avoid  the  storage  of 
large  amounts  of  time-domain  data  [4]. 

Implementation  of  the  FDTD  simulation  using  MPI 

Modem  compilers  on  computers  with  multiple  CPUs  will  parallelize  computer 
programs  using  an  option  called  Open  MP.  Open  MP  distributes  the  code  among 
the  available  CPUs  in  a  computer.  It  will  not  distribute  the  code  among  different 
computers. 

Message  Passing  Interface  (MPI)  is  a  software  package  that  allows  the  programmer  to 
decide  how  the  computation  is  distributed  among  the  CPUs  in  one  machine,  or  among 
the  CPUs  in  several  machines  [5-8].  The  implementation  using  MPI  requires 
considerable  additional  programming  effort.  This  was  done  in  the  hope  of 
surpassing  the  speed  achieved  by  Open  MP,  as  well  as  acquiring  the  ability  to 
distribute  the  program  over  several  computers. 
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MPI  utilizes  the  Domain  Decomposition  Method  (DDM),  a  protocol  that  solves  large 
numerical  boundary  value  problems  by  splitting  the  main  problem  space  into  smaller 
sub-domains.  Each  sub-domain  retains  the  original  qualities  of  the  main  domain  and 
follows  the  same  coordinate  structure  as  the  original  problem  space.  The 
non-overlapping  DDM  requires  some  communication  among  the  sub-domains, 
creating  the  need  for  a  message  passing  interface  like  the  MPI  library’.  Figure  1.1 
shows  an  example  of  a  problem  space  that  is  1003  cells  that  can  be  divided  into  eight 
sections  of  503  sub-domains. 


x 


Figure  1.1.  Domain  decomposition  into  eight  independent  sub-domains. 


The  setup  illustrated  in  Fig.  1.1  is  as  balanced  as  possible  and  creates  a  parallel 
environment  so  each  sub-domain  can  communicate  to  the  same  number  of  neighboring 
sub-domains.  This  minimizes  the  lag  that  can  occur  when  a  core  has  to  wait  for  data 
from  another  to  continue  processing.  The  number  of  sub-domains  depends  on  the 
number  of  sections  on  each  side  of  the  main  domain.  A  sub-domain  that  has  all  six 
neighbors  is  illustrated  in  Fig.  1.2. 


Figure.  1 2.  An  illustration  of  how  one  sub-domain  sees  its  six  surrounding  sub-domains. 


Decomposition  of  the  PML  and  the  lossy  media  is  related  to  the  position  of  each 
sub-domain.  Figure  1.3  illustrates  a  two-dimensional  domain  that  has  been  divided 
into  four  sub-domains,  where  each  contains  different  media  and  PML.  The  process 
takes  in  the  topographic  data  in  the  format  of  the  main  domain  by  each  processing  core 
and  keeps  only  the  data  that  applies  to  its  related  sub-domain.  The  algorithm  uses  the 
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position  data  to  realize  its  correct  place  among  other  sub-domains  and  compares  that  to 
the  main  domain’s  three-dimensional  mesh. 


Figure  1.3.  Two-dimensional  illustration  of  the  PML  and  Media  defined  for  4  sub-domains. 


The  Message  Passing  Interface  Library  makes  communication  at  the  sub-domain 
boundaries  possible.  Communication  between  the  sub-domains  is  facilitated  by  MPI’s 
blocking  receive  operation  (MPI  Recv)  and  MPI’s  non-blocking  send  operation 
( MPI  iSend).  The  blocking  receive  operation  keeps  each  sub-domain  from  continuing 
with  the  calculation  of  the  E  and  H  fields,  until  the  needed  parameters  are  received. 
However,  the  non-blocking  send  operation  sends  the  data  to  the  addressed  sub-domain 
and  continues  with  the  rest  of  the  FDTD  calculation.  Using  a  combination  of 
non-blocking  send  operations  and  blocking  receive  operations  reduces  the  chance  of 
possible  software  hang-ups.  In  this  case  hang-ups  occur  either  when  a  message  is 
needed  that  has  not  been  sent  or  a  sub-domain  is  awaiting  a  successful  sent 
confirmation  that  has  not  been  received  by  the  other  sub-domain.  Figure  1 .4  shows  the 
fields  that  must  be  passed  to  the  neighboring  sub-domains. 


Figure  \  4.  Field  information  transferred  at  the  boundaries  of  sub-domains 
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The  following  sequence  insures  that  the  parallel  FDTD  algorithm  performs  smoothly. 


Step  1)  Send  the  correct  E-fields  to  the  back,  left,  and  bottom  sub-domains. 

Step  2)  After  receiving  the  required  E-fields  from  the  front,  right,  and  top 
sub-domains,  calculate  the  H-fields. 

Step  3)  Send  the  correct  H-fields  to  the  right,  top,  and  front  sub-domains. 

Step  4)  Update  the  D-fields.  The  sub-domains  must  receive  the  required  H-fields  from 
the  back,  left,  and  bottom  sub-domains. 

Step  5)  Calculate  the  E-fields  from  the  D-fields. 

As  a  result  of  the  parallel  FDTD  calculation,  the  electromagnetic  field  data  are  defined 
for  each  sub-domain  and  each  sub-domain  can  be  simulated  separately.  At  the  end  of 
the  parallel  FDTD  process  in  each  sub-domain  the  resulting  E  and  H  fields  are  present. 
The  last  sub-domain  accumulates  all  the  fields  in  one  large  array.  These  E  and  H  fields 
will  look  similar  to  results  from  the  regular  FDTD  method  after  executing  sequentially 
across  only  one  core.  This  is  illustrated  in  Fig.  1 .5. 


Figure  I  5.  A  three-dimensional  view  of  the  E-field  in  sub-domain  1  (left)  and  the  entire 
domain  (fight). 

Once  the  FDTD  code  had  been  programmed  using  MPI,  a  comparison  was  made  with 
a  similar  program  that  only  used  the  Open  MP.  The  program  problem  space 

was (120)  cells.  It  was  distributed  among  27  cores.  A  total  of  30,000  time  steps 

were  needed.  Open  MP  preformed  the  simulation  in  15  minutes,  56  seconds.  The 
MPI  code  performed  the  same  simulation  in  9  minutes,  54  seconds.  Therefore,  MPI 
reduced  the  computation  time  by  35  %  [1,  2]. 

Implementation  of  the  FDTD  code  using  dedicated  hardware. 

Another  approach  to  increasing  the  speed  of  very  large  FDTD  programs  was  the 
implementation  on  dedicated  hardware  using  a  system  from  Acceleware  Corporation. 
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Accelevvare  systems  exploit  the  performance  of  video  graphics  cards  to  obtain 
substantial  speed-up  (Fig.  2.1)  [9]. 


HP  Workstation 
XW9400 


/\ 

\/ 


Software  Developement 
Kit  (SDK)  Version  9.2.0. 


Graphics  Processing  Unit  (GPU) 
NVIDIA  Quadro 
FX  5600  video  card 


Figure  2. 1  Block  diagram  of  the  Acceleware  system 

The  Software  Development  Kit  (SDK)  is  provided  by  Acceleware.  It  is  a  high-level 
programming  language  written  in  C++.  The  FDTD  codes  must  be  rewritten  in  this 
language.  This  language  incorporates  “handles’'  which  represent  different  parts  of 
the  FDTD  simulation.  Table  2.1  is  a  list  of  some  of  the  important  handles  and  their 
functions.  Figure  2.2  is  a  flow  chart  showing  the  steps  for  the  implementation  of  an 
FDTD  program  using  Acceleware. 

Table  2. 1 .  The  Handles  used  in  SDK 


Handle 
Ax_timeexc_t 
Ax_rgnhandle_t 
Ax_mathandle 
Ax  simhandle 


Type 

Time  excitation 
Region 
Material 

Simulation  handle 
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Programming  Flow 


Figure  2.2.  The  programming  How  to  implement  an  FDTD 
program  in  Acceleware 

A  test  was  made  to  insure  that  the  Acceleware  FDTD  was  in  agreement  with  the 
standard  FDTD.  The  simulation  problem  space  and  the  results  are  shown  in  Fig.  2.3. 


monitor 


800  m  Distance  (m) 

(a)  (b) 

Figure  2.3,  (a)  Test  configuration  (b)  Previous  FDTD  (solid  line)  vs.  Acceleware  calculation  (circles). 

Table  2.2  summarizes  the  wall  clock  times  needed  for  an  FDTD  simulation  of  (120) 
cells  over  30,000  time  steps. 

Table  2.2.  The  wall  clock  times  required  for  an  FDTD  simulations 

of(l20)  cells  for  30,000  time  steps.  The  first  three  entries  were  using 

the  8  Quad-core  AMD  Opteron™  Processor  8380.  The  Acceleware 
simulation  was  done  with  Acceleware  libray  version  9.2.0  and  a  NVIDIA 
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Quadro  FX  5600  video  card. 


1  core  63  m  18s 

32  core  (Open  MP)  12m  56s 

27  core  (MPl)  9m  54s 

Acceleware  4m  32s 


Time-domain  Near  to  Far  Field  Transformation 

A  near-to-far  field  transformation  had  previously  been  developed  [10].  This 
transformation  was  developed  to  address  the  problem  of  having  to  model  an 
electromagnetic  source  like  a  ship  at  one  a  relatively  small  resolution  like  10  meters, 
and  yet  having  to  model  a  very  large  problem  space  where  larger  resolution  like  50 
meters  would  be  more  desirable  (Fig.  3.1)  The  previous  formulation  of  this  near-to-far 
field  transformation  stored  the  amplitude  and  phase  of  selected  frequencies  of  each  of 
the  E  fields  on  the  equivalence  surface.  This  amplitude  and  phase  was  used  in 
generating  the  sinusoidal  source  for  the  far  field.  The  disadvantage  of  this  approach 
is  that  a  separate  simulation  is  needed  at  every  frequency  of  interest  in  the  far  field. 
If  instead  the  time  domain  data  at  the  fields  on  the  equivalence  surface  could  be  stored, 
then  information  at  all  frequencies  of  interest  could  be  obtained  in  the  far  field  with 
just  one  simulation.  The  problem  is  that  the  time  domain  data  in  the  near  field  could 
be  several  thousand  points,  and  storing  all  this  data  at  all  the  fields  of  interest  is  not 
practical. 


(a)  It  might  be  desirable  to  model  the  EM  radiation  of  a  ship  over  very  large 
distances. 
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I - T 


(c)  Far  Field  (50  rrf  cells) 


Figure  3.1.  (a)  A  projected  simulation  problem  where  it  would  be  desirable  for  a  source 

of  EM  radiation,  such  as  a  ship,  to  be  modeled  with  relatively  small  cells  of  about  ten 
meters  squared.  However,  problem  spaces  ranging  over  distances  of  three  or  more 
kilometers  would  require  too  many  cells,  (b)  A  near  field  calculation  determines  the 
radiation  from  the  source  at  an  equivalence  surface  surrounding  the  source,  (c)  This 
equivalence  surface  is  then  used  as  the  source  in  an  FDTD  problem  with  much  larger  cells. 

In  order  to  minimize  the  amount  of  data  that  has  to  be  stored  in  the  near  to  far  field 
transformation,  wavelets  are  being  used  to  compress  the  data  [11].  The  technique  of 
compressing  data  in  an  FDTD  simulation  has  been  used  by  this  research  group 
previously  [12,  13].  The  type  of  wavelet  processing  structure  being  used  is 
illustrated  in  Figure  3.2.  The  circles  with  arrows  pointing  downwards  indicate 
“down-sampling,”  i.e.,  every  other  data  point  is  eliminated.  The  circles  with  arrows 
pointing  up  indicate  “up-sampling,”  i.e.,  a  zero  is  added  after  each  data  point.  The 
squares  indicate  convolution  with  the  filter  written  inside.  One  such  group  of  filters 
is  shown  in  Figure  3.3. 
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(a)  Analysis  tree 


Figure  3,2.  The  type  of  structures  being  used  to  achieve  data  compression  and 
reconstruction  using  wavelets.  The  circles  with  arrows  pointing  downwards  indicate 
“down-sampling,”  i  e.,  every  other  data  point  is  eliminated.  The  circles  with  arrows 
pointing  up  indicate  “up-sampling,”  i.e.,  a  zero  is  added  after  each  data  point.  The 
squares  indicate  convolution  with  the  Filter  written  inside.  Each  stage  in  the  analysis  tree 
separates  the  data  into  low  pass  parts  (cOs)  and  high  pass  parts  (els).  Note  that  it  is  only 
the  low  pass  part  that  is  processed  further  The  values  obtained  from  the  analy  sis  can  be 
used  to  reconstruct  the  original  waveform  in  the  synthesis  tree. 
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Figure  3.3.  Four  filters  of  the  type  used  by  the  structure  shown  if  Fig.  3.2.  The 
filters  ho  and  f0  are  low  pass  filters  while  h|  and  fj  are  high  pass  filters. 

As  an  example,  Fig.  3.4  shows  a  waveform  of  6000  time  steps.  An  analysis  tree  of 
eight  levels  is  used  to  produce  the  data  shown.  The  low  pass  components  are  the 
solid  lines  while  the  high  pass  components  are  the  dashed  lines.  Notice  that  the 
dashed  lines  are  virtually  zero  compared  with  the  solid  lines.  Therefore,  all  high 
pass  data  can  be  discarded  and  the  original  waveform  can  be  reconstructed  from  the 
forty-five  points  of  level  eight  as  shown  in  Fig.  3.5. 

Since  we  have  determined  that  only  the  low-pass  parts  of  the  analysis  and  synthesis 
trees  are  needed,  the  implementation  of  the  analysis  or  synthesis  only  requires  about 
ten  lines  of  additional  code  in  the  FDTD  simulations. 
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Figure  3.4.  A  waveform  with  6000  points  (top,  left).  The  subsequent  plots  show  the  outputs  of  the 
various  levels  in  the  analysis  tree.  By  level  8  (lower  right),  only  forty-five  points  are  needed  to 
represent  the  original  signal. 


Time  steps 

Figure.  3.5  The  reconstructed  waveform  (dashed  line)  and  the  original  (solid  line). 

Figure  3.6  is  an  illustration  of  a  simulation  to  show  that  the  near-to-far  field 
transformation  gives  the  same  results  as  a  corresponding  FDTD  program  with  no 
transformation.  The  source  is  a  magnetic  dipole  just  below  the  water  surface.  The 

near  field  program  is  (120)  cells.  Each  cell  is  ten  meters  cubed.  The  far  field 
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simulation  is  only  (60)  cells,  but  the  cells  are  fifty  meters  cubed.  Comparisons 

are  made  at  three  points  which  are  a  lateral  distance  of  500  m  from  the  dipole.  Point 
A  is  100  m  above  the  water  surface,  point  B  is  100  m  below  the  surface,  and  point  C 
is  300  m  below  the  surface.  The  results  are  shown  in  Fig.  3.7. 


The  source  is  a 
dipole,  5  m  below 
the  water  surface 


]QOm 
-100  m 
-300  ni 


500  m 

Figure  3.6.  Diagram  of  the  near  field  problem  space  used  to  shown  that 
the  field  produced  in  the  far  field  are  the  same. 


Figure  3.7.  Results  of  the  test  described  in  Fig.  3  6.  T  he  waveform  on  top  is  source  at  the  dipole 
The  next  three  plots  compare  the  near  and  far  field  Hz  fields  at  points  A,  B,  and  C.  The  solid  lines  are 
from  the  near  field  calculation;  the  dashed  lines  are  from  the  far  field  calculation 
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Appendix  D 


Multiple-layered  Quasi-Electrostatic  (QES)  Development 


Robert  Rebich,  Jeffrey  Young  and  Chris  Wagner 


Department  of  Electrical  and  Computer  Engineering,  University  of  Idaho 


Analytical  Development 


Assume  a  domain  where  the  media  is  composed  of  simple  matter,  in  which  case 


£2 

ii 

Q 

(1) 

B  =  pH 
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Jc  =  crE. 
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Here  D  is  the  electric  displacement  density,  B  is  the  magnetic  flux  density,  Jc  is  the 
electrical  conduction  current  density,  E  is  the  electric  field  intensity  and  H  is  the 
magnetic  field  intensity.  The  permittivity  e  is  a  product  of  the  relative  and  free  space 
permittivity  so  that  c  =  cr e0,  where  f0  =  8.854  x  10-12  F/m.  The  domain  is  absent 
of  all  magnetic  materials  so  that  the  permeability  is  equal  to  that  of  free  space, 
p  =  fjLo,  where  p0  =  x  10-7  H/m.  The  electrical  conductivity  is  represented  by 


The  fields  within  the  domain  are  deemed  quasi-electrostatic  when  the  magnetic 
field  has  little  to  no  time  variation  such  that 


(4) 


As  a  consequence  of  Eqn.  (4),  Faraday’s  law  states  that  the  curl  of  the  electric  field 
is  then  approximately  zero,  in  which  case 


VxE^O. 


(5) 


From  Ampere’s  law, 


V  x  H  = 


(6) 
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where  J  represents  the  combination  of  conduction  and  impressed  current  densities: 


J  =  Jc  +  J'.  (7) 

By  taking  the  divergence  of  Eqn.  (6)  and  knowing  that  the  divergence  of  a  curl  is 
always  zero,  we  find  that 

v'  {w  +  J)  =V-(VxH)=°.  (8) 

Given  Eqn.  (5),  the  electric  field  at  a  given  point  in  space  is  equal  to  the  negative 
gradient  of  the  electric  scalar  potential  V  at  that  point; 

E  =  -  VK  (9) 

For  homogeneous  media,  it  follows  from  Eqns.  (8)  and  (9)  and  from  the  constitutive 
relationships  of  Eqns.  (1)  and  (2)  that 

4vV  +  (tVV  =  VJ',  (10) 

at 

where  V2  is  the  Laplacian  operator.  The  equation  of  continuity  states  that, 

=  <"> 

where  p  is  the  impressed  charge  density,  so  that 

eg*V  +  ,W~g.  ( ,2) 

In  the  frequency  domain,  equation  Eqn.  (12)  is  similar  to 

vV  =  — ML.  (13) 

a  +  jujt 

where  an  e+JU,t  time  factor  is  assumed.  A  special  note  is  made  that  V  and  p  in  Eqn. 

(12)  are  referenced  in  the  time  domain  (i.e.  V  =  V(t),p  —  p(t))  and  V  and  p 

in  Eqn.  (13)  are  referenced  in  the  frequency  domain  (i.e.  V  =  V(u>),p  =  p(u>)). 
Subsequent  analysis  will  be  restricted  to  the  frequency  domain  so  that  no  ensuing 
confusion  should  remain. 

Suppose  we  have  an  interface  of  two  dissimilar  media  according  to  Figure  1 .  It 
then  follows  from  Eqn.  (5)  that 


(Eyi  -  Efl)  X  n  =  0. 


(14) 


The  total  induced  current  within  a  specific  region  is  given  by, 

J*  =  (a  +  jujt)  E. 


(15) 
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Region  B 

Figure  1 :  A  depiction  showing  n  for  a  two-layered  geometry 


Given  Eqn.  (15),  continuity  of  normal  current  states  that 

n  •  =  n  J'u 


(16) 


or, 

(a A  +  Ma)  n  •  Ea  =  (ctb  +  juts)  n  Efl.  (17) 

Once  V  is  determined  from  solving  Eqn.  (13)  in  the  context  of  the  boundary  condi¬ 
tions  of  Eqns.  (14)  and  (17)  we  then  use  Eqn.  (9)  to  determine  the  electric  field  E. 

Let  us  now  consider  a  single  point  charge  of  strength  q  located  at  the  origin  in 
unbounded  media.  The  electric  potential  is  a  solution  to  Eqn.  (13)  such  that 


V  = 


M 

4n(a  +  jut)r 


(18) 


where  _  _ 

r  =  s/x2  +  y2  +  ~2  -  vV2  +  ~2-  (19) 

Eqn.  (18)  is  known  as  the  Green’s  function  solution  of  a  point  charge  at  the  origin 
in  a  lossy  homogeneous  medium.  This  solution  can  be  equally  expressed  in  integral 
form  by  noting  that  =  0,  in  which  case  Eqn.  (13)  is  equivalent  to 


ld_ 
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J^P 

<7  -I-  juje 


(20) 


The  solution  to  Eqn.  (20)  is  a  combination  of  Bessel  and  exponential  functions: 


V  = 


juq 


47T  (cr  +  jut) 


roc 

/  Jo  M 

J  0 


e~x^d\. 


(21) 


Now  if  the  charge  is  located  at  z  =  h,  Eqn.  (21)  may  be  written  as 


1/  = 


juq 


47t  (<t  4-  jut) 


roc 

/  J0(Ap)e-A|z-/,|dA. 

Jo 


(22) 


With  the  potential  determined  for  a  charge  in  a  single  homogeneous  media, 
the  analysis  can  be  further  extended  to  a  three-layered  media  problem  depicted  in 
Figures  2.  Charges  and  observations  in  region  3  will  be  ignored  for  the  remaining 
development. 
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Figure  2:  Point  charge  in  region  1  for  a  three-layer  configuration. 
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Figure  3:  Superposition  of  potentials  for  regions  1  and  2. 
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Source  Charge  in  Region  1 


According  to  Figure  2,  the  charge  is  placed  in  region  1  at  a  distance  h  above  the 
2  =  0  interface.  We  construct  potential  solutions  for  each  region  as  follows.  For 
observation  locations  in  region  1,  where  2  >  0,  the  total  potential  at  any  given  point 
is  the  superposition  of  two  individual  potentials  as  shown  in  Figure  3a.  The  first  po¬ 
tential  represents  the  direct  path  from  the  charge  to  observer  and  is  of  unity  strength. 
The  second  potential  represents  the  reflected  path  from  the  2  =  0  boundary  and  is 
of  strength  R.  The  total  potential  at  any  location  in  region  1  is  given  by, 

Vi  =  J™  [e-A|z-/l1  +  Re-X^)  J0(\p)dX  (23) 


where, 


Y\  —  0\  +  juJti, 


(24) 


which  is  the  admittivity  of  region  1 .  For  observation  locations  in  region  2,  where 
2  <  0  and  2  >  -d,  the  total  potential  at  any  given  point  is  the  superposition  of 
two  individual  potentials  as  shown  in  Figure  3b.  The  first  potential  represents  the 
downward  traveling  path  caused  by  the  transmission  of  the  direct  path  through  the 
2  =  0  boundary  and  is  of  strength  B.  The  second  potential  represents  the  upward 
traveling  path  caused  by  the  downward  path  reflection  at  the  2  =  —  d  boundary  and 
is  of  strength  A.  The  total  potential  at  any  location  in  region  2  is  given  by, 


r  00 

Vi=A^vr  /  [Ae~X{z+h)  +  J0{. Xp)dX.  (25) 

<^7r  *1  Jo 


For  observation  locations  in  region  3,  where  2  <  -d,  the  potential  at  any  given 
point  represents  the  transmitted  path  of  the  downward  traveling  path  in  region  2  as 
it  encounters  the  z  =  -d  boundary  and  is  of  strength  T.  This  situation  is  shown  in 
Figure  3b.  The  total  potential  at  any  location  in  region  3  is  given  by, 

•  r  co 

V3  =  XxJo  TeX{i~h)j°MdX-  (26) 


Now  that  the  potential  solutions  are  formulated,  the  unknowns  coefficients  R,  A , 
B ,  T  are  found  by  applying  the  boundary  conditions  of  Eqns.  (14)  and  (17)  in 
the  context  of  Eqn.  (9).  The  boundary  conditions  must  be  applied  at  both  region 
interfaces,  i.e.  the  2  =  0  and  2  =  -d  boundaries.  Hence, 


V, 


2=0 


H=o 
V3 


z=— (i  ’ 


(27) 


and 


1  dz 


2=0 


2—0 


(28) 
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Applying  the  above  boundary  conditions  to  the  potentials  of  Eqns.  (23)  -  (26),  we 
obtain  a  set  of  four  equations  with  four  unknowns  given  by, 


1  +  R  =  A  +  B 

(29) 

AeXd  +  Be~xd  =  Te~Xd 

(30) 

(31) 

Be~Xd  -  AcXd  =  ~Te~Xd. 

*2 

(32) 

After  many  algebraic  steps,  we  find  that. 

R23e~2Xd  -  R2l 

1  -R23R2le~™ 

(33) 

,  Rail  -  R2i)e~2Xd 

1  -  RnRne-™ 

(34) 

o  (1  -  Bj,) 

1  -  RizRne-™ 

(35) 

(1  —  R2  i  )  ( 1  +  R23) 

1  -  R23R2lc~*Xd  ’ 

(36) 

where  the  interfacial  reflection-like  coefficients  are  given  by. 

R  ->Wl 

R2'  ~  Y,  +  Y2 

(37) 

R  -  -  y3 

fl“-y2  +  n' 

(38) 

The  potential  integrals  for  charge  in  region  1  as  stated  by  Eqns.  (23)  -  (26)  are 
now  fully  specified  and  can  be  computed  numerically,  or  by  image  summations,  as 

described  next. 


Source  Charge  in  Region  1,  Observation  in  Region  1 

In  this  section  we  will  take  the  rather  complex  potential  integral  of  Eqn.  (23)  and 
simplify  it  into  the  form  of  an  infinite  summation.  This  is  necessary  because  the 
integral  will  eventually  be  solved  in  a  numerical  fashion  and  the  form  of  Eqn.  (23) 
can  be  difficult  to  integrate  numerically.  The  following  procedure  will  make  the 
numerical  solution  quick,  efficient  and  very  robust. 

The  potential  expression  of  Eqn.  (23)  is  first  separated  into  two  parts: 

roo  r  oo 

V\  =  I<\  /  e-xlz-h'j0(\p)d\  + K\  Re-X{!+h)MXp)d\,  (39) 
Jo  Jo 
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where  K\  is  the  normalized  point  charge  strength  in  region  1  defined  by, 

juq 


K  ,= 


47ryi 


(40) 


The  first  term  of  Eqn.  (39)  can  be  equally  represented  in  closed-form  by  comparing 
Eqns.  (18)  and  (21),  in  which  case, 


K , 


TOO 

/  e-x'z~h[J0 
Jo 


(Xp)dX  =  — . 

r 


In  the  present  context,  r  is  defined  by. 


with. 


r  =  \//72  +  (2  —  /i)2 


P=  \/^2  +y2- 


(41) 


(42) 


(43) 


Eqn.  (41)  is  commonly  referred  to  as  the  Weber  integral  [2].  Equation  (39)  is 
equally  stated  as, 

K  f°° 

V1  =  —  +  Ki  /  Re~X{z+h)  J0{Xp)dX.  (44) 

r  J  0 

We  next  insert  7?  from  Eqn.  (33)  into  our  expression  to  obtain, 

V'  =  T  +  K'  [  ***>*»)*■  <45> 

This  integral  can  be  expressed  as  an  infinite  summation  by  the  following  Taylor 
series  expansion, 

1  00 

- - =  xn  for  |x|  <  1.  (46) 

1  —  x  ' 

n=0 

This  allows  us  to  rewrite  Eqn.  (45)  as, 

r/  00  pOO 

Vl  =  -l  +  h\T  R$ 3^,  /  (R23e~2Xd  -  Rtl)  e-x^+h+2nd^Jo(XP)dX.  (47) 

r  n=0 


A  close  inspection  of  the  previous  integral  shows  that  it  has  a  closed-form  solution 

[l]; 

rx  .  1 

(48) 


r°°  1 

/  e-^„(Ap)rfA=  - - 

7o  V  a2  +  p2 


From  this  integral  identity,  it  follows  that  Eqn.  (47)  is  equivalent  to 


K 


V,  =  —  +  A-,  £  am 


n=0 


7?23 


R 


21 


where, 


Dx  —  z  +  h  +  2d(n  +  1) 


(49) 


(50) 
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Z?2  —  z  ~\~  h  2dn. 


(51) 


The  infinite  summation  in  Eqn.  (49)  allows  us  to  view  the  potential  at  some  obser¬ 
vation  location  as  an  infinite  summation  of  equivalent  charges  at  different  locations 
with  different  strengths  and  phases.  Figure  4  shows  the  corresponding  image  lo¬ 
cations  and  strengths  represented  by  Eqn.  (49).  The  fully  filled  black  circle  with 
strength  q  represents  the  original  charge  and  is  located  at  height  h.  Images  ao  and  b0 
represent  the  first  and  second  terms  from  the  summation  of  Eqn.  (49)  when  n  =  0 
and  are  presented  as  gray  shaded  circles.  The  image  depths  are  shown  relative  to 
the  z  =  0  interface.  The  remaining  image  terms,  a„  and  bn  represent  the  infinite 
number  of  concurring  images  and  are  presented  as  light  gray  circles  with  dashed 
outlines.  According  to  Eqn.  (49),  we  see  that  the  strengths  for  the  corresponding 
weighted  images  are  as  follows, 


an  —  A'i/?23 


(52) 


bo  —  A'i/?2i 

(53) 

On  =  AlA^X 

(54) 

bn  =  KiR^RZt'- 

(55) 

It  is  insightful  to  note  that  the  image  charges  are  proportional  to  the  original  charge 
q.  For  example,  we  can  take  Eqn.  (52)  and  substitute  K\  from  Eqn.  (40)  to  obtain 


Oq 


4t tYi 


(56) 


where,  W  can  be  defined  as  a  weighting  term. 


W 

vv  a0 


juR  23 
47TV,  ' 


(57) 


It  is  now  obvious  that  each  image  term  is  of  strength  q  multiplied  by  a  complex 
weighting  term  W .  This  suggests  that  the  image  charges  are  out  of  phase  with 
the  original  charge  and  scaled  appropriately.  It  is  also  important  to  note  that  as  n 
increases,  the  weighting  terms  become  increasing  small  and  the  distances  relative  to 
the  observation  point  become  increasingly  large.  This  means  that  the  summation  of 
the  images  converges  very  rapidly.  No  more  than  a  few  terms  from  the  summation 
are  needed  for  accurate  results  when  dealing  with  any  combination  of  charge  and 
observation  locations. 

The  electric  field  vector  can  be  determined  by  taking  the  gradient  of  the  potential 
Eqn.  (49), 


E  =  -VV  =  - 


/  dV  „  dV  .  0V„  \ 

(  &  a' +  V" +  av  • 


(58) 


The  electric  field  components  are  hence  given  by, 


Ex=I^x  +  I<1xY,R'kRZx 


n=0 


R, 


23 


^21 


n/^F+7  y/Dl  +  7  J 


(59) 
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E,  =  ^z  +  K,Y.^, 


n= 0 
oc 


A 


23 


/?■ 


2t 


n=0 


.y/D  i +  7 

^23^1 


yor+7  j 


^21  D'Z 


(60) 

(61) 


Lv^-k?  n/^I  +  7  J 

Source  Charge  in  Region  1,  Observation  in  Region  2 


The  aforementioned  procedure  may  also  be  applied  to  Eqn.  (25)  to  change  the 
complex  integral  into  an  infinite  summation.  The  detailed  steps  will  not  be  shown 
for  this  equation  due  to  the  similarity  with  the  previous  process.  The  final  form  of 
the  summation  equation  is  given  by. 


V2  =  Ky{\  —  R21)  ^2  ^23-^21 


n=0 


R 


23 


+ 


1 


where, 


£>3  =  h  —  z  +  2  dn. 


(62) 


(63) 


Figure  5  shows  the  corresponding  image  locations  and  strengths  represented  by 
Eqn.  (62).  The  strengths  for  the  corresponding  weighted  images  are  as  follows. 


Oo  —  AT(1  —  R21 ) 

bo  =  AT  (1  —  A2OA23 

01  =  AT(1  —  A2OA23A21 

bi  —  I<i(l  —  A2i)A|3A2i 
an  =  A  i(l  —  A21 )  R^iR^zi 

bn  =  AT(  1  -  A2i)A2n3+X- 

The  electric  field  components  are  hence  given  by, 

OC 

Ex  =  ATx(1-A2I)^A53A2T 

n=0 

oc 

hy  =  K\y(  1  —  R 21)  'y  "  R*23^2\ 
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v/Bf  +  Z 

^0 

to 

00 

b 

1 

1 

Dz 

v/^f+7 


(64) 
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Figure  5:  A  depiction  of  the  image  principle  per  Eqn.  (62). 
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Source  Strength  K\ 


In  this  section  the  point  charge  strength  Kx  will  be  specified  in  terms  of  a  two  elec¬ 
trode  problem.  We  will  begin  with  the  definition  of  K\  from  Eqn.  (40).  We  will 
assume  the  existence  of  two  spherical  electrodes  of  charge  +  q  and  -q,  of  radius  a 
and  of  separation  L.  We  will  define  the  configuration  of  two  oppositely  charged 
electrodes  as  a  source.  If  we  assume  the  charge  distribution  is  uniform  across  the 
entire  surface  of  both  electrodes  then  Eqn.  (49)  is  valid  for  the  following  proce¬ 
dure.  Given  the  charge  separation,  a  potential  Vo  is  assumed  to  exist  between  the 
electrodes.  Let  P i  be  any  point  on  the  negative  charged  sphere  and  P2  be  any  point 
on  the  positive  charged  sphere.  Then  by  Eqn.  (49)  for  source  location  in  region  1 
and  from  superposition,  we  have, 

P2 

(73) 

Pi 

We  will  define  the  geometrical  factor  G\  to  be  the  following: 


Vo  =  Kx 


1  00 

- + E 


n=0 
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___R 23_ 

VD\  +  p2 


R 


21 


yM+72 


Pi 
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(74) 


The  source  strength  is  now  determined  by  supplying  the  potential  between  the  two 
electrodes  and  calculating  the  geometrical  factor; 


Kx 


Vb 

G\ 


(75) 


Source  Charge  in  Region  2 


In  Figure  6  the  charge  is  located  in  region  2  at  a  distance  —  h  below  the  2  =  0 
interface.  We  construct  potential  solutions  according  to  Eqn.  (22)  for  each  region 
as  follows.  For  observation  locations  in  region  1,  where  2  >  0,  the  total  potential 
at  any  given  observation  point  is  equal  to  the  upward  transmitted  path  through  the 
2  =  0  boundary,  which  is  of  strength  Tx  and  is  shown  in  Figure  7a.  The  total 
potential  is  given  by, 


roo 

V'  =  TKT  Txe-XzJ0(XP)dX.  (76) 

47tY2  Jo 

For  observation  locations  in  region  2,  where  2  <  0  and  2  >  —  d,  the  total  potential  at 
any  given  observation  point  is  equal  to  the  superposition  of  three  different  potentials 
and  is  shown  in  Figure  7b.  The  first  is  the  direct  path  from  charge  to  observer  and  is 
of  unity  strength.  The  second  is  the  upward  traveling  path  caused  by  the  reflection 
from  the  direct  charge  from  the  2  =  —d  boundary  and  is  of  strength  A.  The  last 
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Figure  6:  Point  charge  in  region  2  for  a  three-layer  configuration. 
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Figure  7:  Superposition  of  potentials  for  regions  1  and  2. 
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is  the  downward  traveling  path  caused  by  the  reflection  of  the  direct  charge  off  the 
2  =  0  boundary  and  is  of  strength  D.  The  total  potential  is  given  by, 


V2  =  p^r  1  +  Ae-Az  +  Bcx*} J0(Xp)dX. 

471*2  JQ 

(77) 

For  observation  locations  in  region  2,  where  2  <  —  d,  the  total  potential 
the  direct  charge  transmitted  path  through  the  z  —  -d.  boundary  and  is 
T2,  and  is  shown  in  Figure  7b.  The  total  potential  is  given  by, 

is  equal  to 
of  strength 

V3  =  TV  [*  T2ex*J0(\p)d\ 

47rV2  Jo 

(78) 

The  same 
boundary 

boundary  conditions  of  Eqns.  (27)  and  (28)  still  apply.  Applying  these 
conditions,  we  obtain  a  set  of  four  equations  with  four  unknowns: 

T\  =  exh  +  A  +  B 

(79) 

T2e~xd  =  e“A(/l+d)  +  AeXd  +  Be~Xd 

(80) 

Y\T\  =  Y2{cxh  +  A  -  B) 

(81) 

Y2(AeXd  -  Be.-Xd  -  e~X{h+d))  =  -Y3T2e-Xd. 

(82) 

Solving  for  the  unknown  coefficients  of  interest,  we  obtain  the  following  equations: 


(R2l  +  l){exh  +  R23e-X^h+^) 

1  -  R23R2 \e~2Xd 

(83) 

R23(l  +  R2le2Xh)e-x(h+2V 

1  -  R23R2le~^ 

(84) 

R2\exh  +  R23R2le-*h+2V 

1  -  7?237?21e_2Ad 

(85) 

Source  Charge  in  Region  2,  Observation  in  Region  1 

The  image  solution  for  the  configuration  of  Figure  6  with  observation  in  region  1  as 
stated  by  Eqn.  (76)  is  given  as, 


V\  =  A2(l  +  R2l)  ^  A53/?4'i 


71—0 


R 


23 


+ 


\fW+?  y/Dl+7 


where, 


Da  =  z  —  h  +  2  fin 


and  h'2  is  the  point  source  strength  in  region  2  defined  by, 


f<7  = 


4  ttY2 


(86) 

(87) 


(88) 
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Figure  8:  A  depiction  of  the  image  principle  per  Eqn.  (86). 
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Figure  8  shows  the  corresponding  image  locations  and  strengths  represented  by 
Eqn.  (86).  The  strengths  for  the  corresponding  weighted  images  are  as  follows, 


«0  —  ^2(1  +  7?2l) 

bo  =  A  2^23(1  +  7?2i) 

fli  =  /^2i?23(l  +  7?2l)7?21 
61  =  +  7?2l)/?21 

an  =  K2R&(1  +  7?2i)7^1 
()n  =  /\27?^3+,(l  +  /?2i)7^1. 
The  electric  field  components  are  hence  given  by, 

OC 

Ex  =  Klx(l  +  Rn)'jTt  7^3^21 


n=0 


Ev  =  Kly(l  +  R2i)'£iRbRZ\ 

n= 0 
oc 

EZ  =  KX{  l  + 


n=0 


7?23 

1 

Vdi +  p2 

y/D24+ 

%  3  + 

1 

Lv^fT? 

■JD\  +  / 

^1^23 

n4 

.v^f+7'  x/^f+7. 


(89) 

(90) 

(91) 

(92) 

(93) 

(94) 


(95) 

(96) 

(97) 


Source  Charge  in  Region  2,  Observation  in  Region  2 


The  image  solution  for  the  configuration  of  Figure  6  with  observation  in  region  1  as 
stated  by  Eqn.  (77)  is  given  as, 


where. 


K  00 

-i  +  Ki  £ 

'  ti 

7?23 

^23^*21 

WD\  +  (? 

V»25  +  p2 

j  ^23^21 

7(21 

(98) 

■JD'i  +  p2 

x/T^  +  P2] 

D5  =  z  —  h  +  2 d(n  4-1), 

(99) 

Da  =  h  —  z 

+  2  d(n  +  1), 

(100) 

Di  —  2nd  —  z  —  h. 

001) 
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h  +  2nd  +  2d 
2nd  -  h 


-(h  +  2nd  +  2d) 
h  -  2nd  -  2d 


Figure  9:  A  depiction  of  the  image  principle  per  Eqn.  (98). 
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Figure  9  shows  the  corresponding  image  locations  and  strengths  represented  by 
Eqn.  (98).  The  strengths  for  the  corresponding  weighted  images  are  as  follows, 


On  =  7^2^23^21 

bo  =  7\27?237?21 
Co  =  7\27?2i 

/o  =  7^2  7?23 

«n  =  7^27?2n3+17?2n1+1 
6n  =  K2K£xFC£x 
Cn  =  7^7?2n37?2n1+1 
fn  =  7i27?2n3+1772ni- 
The  electric  field  components  are  hence  given  by, 


K- 


Ex  =  -ix  +  K2xJ2^3Rn2i 


n=0 


ll 


23 
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7?237?2i 

Lx/TTNV'5  v/^I+73 
7?237?21  ,  7{21 
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y/Dfr^  y/W+rl 


K. 


Ey  =  -ly  +  K2yJ2R^ 


n= 0 
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23 


y^F+7 
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7?  7? 


23^*21 


+ 


R23R2I 
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y/D\  +  P2 
^21 


\/  7>I  +  P2  v^f  +  P2 


7T 


^  =  -y*  +  7^2  53  7?-3i?” 


n=0 


79i7?: 


1  -*  ^23 


Ly^fTp-2 

Dg/?23/?21 
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D5R23R21 

y/Dl  +  P2* 

D7R. 


21 


VW+7  J 


Source  Strength  I<2 


(102) 

(103) 

(104) 

(105) 

(106) 

(107) 

(108) 
(109) 


(HO) 


(HI) 


(112) 


For  source  location  in  region  2  and  by  Eqn.  (98)  we  have: 


l/o  =  7^ 


1  00 

-  +  E 


n=0 


7?23  7?237?21  7?23  7?21  7?21 

D2+p2  D52  +  P2  7)2+p24  DHvJJ^  • 

(1*13) 
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We  will  define  the  geometrical  factor  G2  for  the  two  electrode  scenario  to  be  the 
following: 


1 

C2  =  -  +  Y. 

n=0 


#23  #23#21  #23#21 

D1  +  P2  Dl  +  p>  Dl  +  p2  " 


#21 

Dj  +  p2 


Pi 


(114) 


We  now  can  determine  the  source  strength  by  supplying  the  potential  between  the 
two  electrodes  and  calculating  the  geometrical  factor; 


I<2  = 


(115) 


Results 


In  this  section  we  will  use  the  closed  form  Sommerfeld  solution  (denoted  as  WSU, 
per  Appendix  B)  to  compare  and  validate  the  QES  solution.  The  Sommerfeld  so¬ 
lution  has  been  fully  validated  and  is  regarded  as  the  exact  full-wave  solution.  For 
the  first  scenario  we  place  the  source  in  region  2  along  the  x-axis  at  a  depth  of  10m 
with  a  normalized  strength  of  1  A-m.  Regions  1,  2  and  3  will  be  that  of  air,  water 
and  mud  respectively,  as  specified  in  Table  1.  The  electric  held  magnitude  is  ob¬ 
served  along  the  horizontal  x-axis  at  a  depth  of  50m.  The  WSU  solution  relies  on  a 
point  source,  so  for  an  equivalent  scenario,  the  QES  source  will  be  made  small.  The 
QES  source  electrodes  have  a  radius  of  0.2m  and  are  separated  by  1  m.  The  results 
for  this  scenario  are  shown  in  Figure  10.  The  QES  solution  has  little  change  with 
frequency  while  the  WSU  solution  is  frequency  dependent.  For  this  scenario  the 
QES  solution  falls  within  50  percent  of  the  WSU  solution  for  frequncies  up  to  1000 
Hz  and  is  very  accurate  at  lower  frequencies.  For  the  next  scenario,  the  same 


Table  1 :  Material  Properties 


Region  1 

Region  2 

Region  3 

Cr  1 

81 

1 

a  0 

0.018 

0.012 

source  location  and  strength  will  be  used.  Region  properties  will  be  consistent  with 
Table  1 .  We  observe  the  electric  field  magnitude  along  the  vertical  2-axis  at  a  radial 
distance  of  100m.  The  results  are  shown  in  Figure  1 1 .  The  results  are  similar  to  that 
of  the  first  scenario  except  more  accurate.  Consider  the  scenario  when  the  source 
is  moved  further  in  the  water  to  a  depth  of  50m  and  the  observation  is  moved  out  a 
distance  of  300m.  The  data  for  this  scenario  is  shown  in  Figure  12.  We  can  see  that 
the  QES  solution  is  highly  accurate  when  the  frequency  is  low.  For  this  scenario, 
when  the  frequency  is  increased  upwards  to  1000  Hz,  the  QES  solution  over  pre¬ 
dicts  the  magnitude  of  the  electric  held  by  up  to  150  percent  but  is  still  precise  at 
low  freqeuncies.  Clearly  the  region  of  validity  for  the  QES  solution  is  determined 
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Figure  10:  HED  in  water  -  Horizontal  sweep 

by  the  relative  distance  between  source  and  observation,  and  frequency.  In  the  near 
field,  the  QES  solution  predicts  to  a  high  level  of  accurately  while  in  the  far  field, 
the  QES  solution  becomes  increasingly  inaccurate.  Also  we  notice  that  the  QES 
solution  error  is  larger  in  water  than  in  air.  Further  work  is  needed  to  ascertain  the 
QES  region  of  validity. 
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E-Field  Magnitude  (V/m)  E-Field  Magnitude  (V/m) 


Vertical  distance  (m) 

Figure  1 1:  HED  in  water  -  Vertical  sweep  100m  away 


Figure  12:  HED  in  water  -  Vertical  sweep  300m  away 
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A  Perfectly  Matched  Layer  for  Lossy  Media 
at  Extremely  Low  Frequencies 
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Abstract — The  perfectly  matched  layer  (PML)  has  proven  to  be 
an  effective  means  of  absorbing  outgoing  waves  for  finite-differ¬ 
ence  time-domain  (FDTD)  simulations.  This  letter  describes  the  de¬ 
velopment  of  a  PML  specifically  for  underwater  simulations  at  low 
frequencies.  This  is  a  significant  development  for  this  project  that 
involves  simulations  of  electromagnetic  signals  for  long  distances 
under  water. 

Index  Terms — finite-difference  time-domain  (FDTD)  methods, 
perfectly  matched  layer  (PML). 


I  Introduction 

THE  single  largest  threat  to  surface  warships  is  mines. 

These  mines  are  often  detonated  by  the  electromagnetic 
signature  of  a  surface  ship  [1].  For  this  reason,  it  is  desirable  to 
have  simulation  methods  to  study  the  propagation  of  extremely 
low  frequency  (ELF)  electromagnetic  waves  under  water.  The 
finite-difference  time-domain  (FDTD)  method  [2],  [3]  is  one  of 
the  most  widely  used  methods  in  electromagnetic  simulation 
and  has  recently  been  adapted  for  ELFs  under  water  [4],  [5]. 
In  FDTD  simulations,  it  is  necessary  to  have  an  absorbing 
boundary  condition  (ABC)  to  truncate  the  problem  space  and 
absorb  outgoing  waves.  One  of  the  most  widely  used  and 
versatile  ABCs  is  the  perfectly  matched  layer  (PML)  [6],  [7]. 
There  has  been  some  activity  in  the  development  of  PMLs  that 
are  effective  in  low  frequency  or  dispersive  media  [8]— [  11]- 
In  this  letter,  we  describe  the  development  of  a  PML  that  is 
specifically  suited  for  very  lossy  media  at  ELFs. 

II.  Implementation  of  the  PML 

A.  Be  render's  PML  in  Free  Space 

Berenger  [6]  assumed  that  any  plane  wave  propagating  in  the 
direction  d  near  the  PML  could  be  broken  up  into  the  part  trav¬ 
eling  perpendicular  to  the  PML  d±  and  the  part  traveling  parallel 
d||  (Fig.  1).  The  two  conditions  for  the  PML  are  the  following: 
1)  It  must  have  the  same  impedance  as  free  space  and  not 
present  a  loss  to  the  wave  traveling  parallel  to  the  interface. 


Free  space 

PML  Medium 

1 

k 

d'j 

A 

Fig.  1.  The  PML  is  implemented  by  assuming  any  propagating  wave  can  be 
broken  up  into  a  part  that  is  perpendicular  to  the  PML  interface  and  a  part  that 
is  parallel  to  it. 


2)  It  must  increase  the  artificial  electric  and  magnetic  conduc¬ 
tivities  such  that  the  impedance  still  matches  that  of  the  free 
space. 

Both  of  these  conditions  are  met  by  increasing  the  electric  and 
magnetic  conductivities  in  the  PML  such  that 


7/  = 


1  + 

ju>n0  J 

1  +  4^- ) 


(I) 


Note  that  this  impedance  is  a  real  number. 

Berenger  implemented  (1)  into  the  FDTD  formulation  by  a 
split-step  formulation  that  broke  each  electric  and  magnetic  field 
into  two  components.  Most  applications  assume  that  the  back¬ 
ground  medium  in  the  main  problem  space  is  free  space. 


B  The  PML  in  a  Lossy  Medium  at  ELFs 

When  the  background  medium  is  lake  water  and  the  frequen¬ 
cies  are  in  the  ELF  region,  the  situation  is  different.  Lake  water 
has  a  dielectric  constant  of  80  and  a  conductivity  of  0.018  S/m 
[  1 2].  For  lake  water  at  1  kHz,  the  complex  dielectric  constant  is 

*  _  0.018 

—  £w  H — : -  —  80  H - - o  \  ~7 - ryr 

jueo  j  (2 7T  x  103)  (8.85  x  HT12) 

=  80  —  j3.24  x  105  S  — j3.24  x  105. 
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Therefore,  the  impedance  is 


This  impedance  can  be  written  in  polar  coordinates  as 

Vw  =  F/u  l  ^45°. 
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Notice  that  because  the  loss  term  dominates,  the  impedance  is 
at  45°.  The  impedance  of  the  PML  material  must  remain  at  this 
value,  but  at  the  same  time  increase  the  loss  further  as  it  goes  per¬ 
pendicular  into  the  PML.  This  can  be  accomplished  by  adding 
a  factor  s  to  the  conductivity  and  the  permeability 


s'lL"  ns 

Vpml  —  .  /  9.gm  •  w) 

V 

This  addition  of  the  a  term  causes  the  PML  medium  to  absorb 
outgoing  waves  faster  than  the  water  medium,  but  also  avoids 
reflections  from  the  PML  medium.  The  factor  s  equals  one  in 
the  background  medium,  but  increases  as  it  goes  into  the  PML. 

This  implementation  may  bare  a  superficial  representation 
to  the  “stretched  coordinates’*  proposed  by  Chew  and  Weedon 
[13].  However,  the  s  in  (3)  is  a  real  number  as  opposed  to  the 
complex  numbers  used  in  the  stretched  coordinates.  There  have 
been  other  methods  proposed  for  the  PML  in  lossy  media  [10], 
[11]  where  the  conductivity  is  large  enough  that  it  plays  a  sub¬ 
stantial  role  in  the  complex  dielectric  constant.  However,  the 
impedance  in  (3)  is  for  the  situation  when  the  imaginary  part  of 
the  dielectric  constant  dominates  completely. 


where  Ax  and  At  are  the  cell  size  and  time-step,  respectively. 
The  implementation  of  the  PML  in  the  II z  field  in  the  y-direc- 
tion  is  relatively  straightforward.  The  term  s  is  added  only  to  the 
y  differential 


=#r 1/2  (*.;>*) 

+  <£L)  t£i‘+1  (*'>*  *>  -  ^l+1  (i  +  i.i. *)] 

(6) 

Instead  of  changing  the  entire  term  containing  s  above,  we  have 
found  it  expedient  to  include  a  one-dimensional  array 

oy  (J)  =  ^  (7) 

so  (6)  becomes 


C.  Implementation  Into  FDTD 

We  begin  with  the  following  formulation  of  the  Maxwell 
equations  for  a  lossy  media: 


0E  ri  „  rr 
^o—*  -  oE  =  V  x  // 
at 

dll  „  „ 


We  will  restrict  the  discussion  to  the  implementation  of  Ex  and 
II z  propagating  in  the  y-direction  perpendicular  to  the  PML 


^u'^0 ' 


0EX 


Ot 


 ohz 

DHy 

(4a) 

+  awEx 

0y 

Oz 
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 0EX 

OEy 

(4b) 

110  dt 

0y 

Ox 

The  FDTD  formulation  leads  to  the  following  coupled  equa¬ 
tions: 


r*n+l 


o 


:E?(i,j,k)  + 


(^u>^o  T  w  ■  At^  T 

II'zl+1/2  (i,j,  k)-Hl,+1/2(i,j-l,k) 

-h:;+1/2  (ij.k)  +  HTyi+1/2  (ij,k- 1) 

EZ+'iiJ  +  hk) -££+,(i, j,k)' 

'  E«+l(i,j,k)-E?+1(i  +  l,j,k) 


At) 

(5a) 


(5b) 


H'zl+1/2  ( i,j,k ) 

=  II’r1/2(iJ,k) 

At 

+  9V  (j)  7 — 

(PoAx) 

•  [£5+1(iJ  +  li*)--B5+1(i,jf*)].  (8) 

Adding  the  PML  to  the  calculation  of  the  Ex  field  requires  that 
the  calculation  be  split  into  two  equations  for  propagation  in  the 
y  and  z  directions.  The  $  term  is  only  added  to  the  conductivity 
in  the  y-direction 


E^(i,j,k) 

=  - rr T£”y  (i,h  k)  +  -r 

•  [//;l+1/2  (i,j,  k)  -  2  (i,  j  -  1,  *)]  (9a) 

EXl(i,j,k) 

-  pn  ( ■  ■  .  x  ,  (S?) 

+  *  ’  (eu  s0  +  ■  AO 

•  [//'i+1/2(z,j,  k  - 1  k)} .  (9b) 

Once  again,  it  is  expedient  to  implement  this  by  an  additional 
one-dimensional  array.  Equation  (9a)  becomes 


E”P  (*>  k)  =  jy  {j)  •  ca  •  Exy  (i,j,  k) 

+fvU )  ■  cb  •  [llll+1/2  ( i.j,k )  -  //"+V2  (i,j  -  1,*)]  (10) 
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80  cells 


Monitor  line 


Ez  Dipole 
source 


1 5  celh 


x  10 


80  cells 

(a) 


(b) 


Fig.  2  (a)  A  dipole  source  is  located  in  the  SO3  problem  space.  Once  steady 

stale  has  been  reached,  the  amplitude  is  determined  at  a  transverse  line  15  cells 
from  the  source.  The  cells  arc  25  m3 .  (b)  The  amplitude  at  the  monitor  line  after 
4000  time-steps. 


where 


CXI  — 

£w£q 

(11a) 

(^lU^O  T  & w  *  A) 

cb  = 

m 

(Hb) 

(£w£q  +  *  A£) 

fy  (j)  = 

(11c) 

( 1  _L_  .  At  V 

In  summary,  the  PML  is  implemented  in  the  y-direction  by 
the  one-dimensional  arrays 


gy  (j)  =  - 


(l  +  aw  ■ 

fy  (j)  =  7^ - ~~t- 


(7*) 

(lie*) 


It  has  been  found  empirically  that  an  effective  formula  for  the 
s  factor  as  it  goes  into  the  PML  is 


a  =  0.2  •  (j  -  >,igo) 


(12) 


where  je<ige  is  the  beginning  of  the  PML.  This  formulation  also 
prevents  the  largest  stretched  cells  from  exceeding  the  skin 
depth. 


III.  Results 

In  this  section,  we  illustrate  the  effectiveness  of  the  lossy 
medium  ELF  PML.  We  will  start  with  the  problem  space  illus¬ 
trated  in  Fig.  2(a),  which  is  80  cells  cubed.  Each  cell  is  25  m 


Fig.  3.  (a)  The  problem  space  is  truncated  to  1 0  cells  to  the  nght  of  the  source. 
A  four-cell  lossy  PML  has  been  added  to  each  boundary  (b)  The  solid  line  is 
the  amplitude  for  the  simulation  in  (a),  while  the  dashed  line  is  from  the  80  cell 
monitor  line  of  Fig.  2,  (c)  The  same  simulation  with  no  PML 


cubed.  The  size  of  80 5  was  needed  so  that  boundary  plays  no 
role  in  the  results  of  the  simulation.  The  source  is  a  single-cell 
electric  dipole.  After  4000  time-steps,  the  amplitude  is  calcu¬ 
lated  via  the  method  of  two  equations,  two  unknowns  (2E2U) 
[14]  at  a  monitor  line  15  cells  from  the  dipole,  as  shown  in 
Fig.  2(b). 

The  simulation  is  then  repeated  for  the  truncated  problem 
space  shown  in  Fig.  3(a),  where  a  four-cell  PML,  has  been  added. 
In  this  simulation,  the  right  wall  has  been  moved  in  to  within  10 
cells  of  the  source.  The  results  are  plotted  in  Fig.  3(b)  (solid  line) 
along  with  the  results  of  the  previous  simulation  (dashed  line). 
For  comparison,  Fig.  3(c)  is  the  same  simulation  with  no  PML 
on  the  truncated  wall.  Clearly,  substantial  errors  appear  when 
the  PML  is  not  present. 

In  one  final  simulation,  the  problem  space  is  reduced  to 
60  x  20  x  20  cells,  as  illustrated  in  Fig.  4(a).  The  results  arc 
shown  in  Fig.  4(b),  where  the  results  of  the  original  simulation 
of  Fig.  2  are  presented  for  comparison.  The  amplitudes  on  the 
monitor  line  within  five  cells  of  the  center  are  identical.  Once 
again,  the  results  without  the  PML  are  shown  in  Fig.  4(c), 
demonstrating  the  expected  error  when  no  PML,  is  present. 
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Fig.  4.  (a)  A  simulation  similar  to  Fig.  2  but  with  the  pmblem  space  truncated 
to  60  x  20  x  20  cells,  (b)  The  solid  line  is  the  amplitude  for  the  smaller  problem 
space,  while  the  dashed  line  is  the  larger  problem  space  of  Fig.  2.  (c)  The  same 
comparison  when  a  60  x  20  x  20  problem  space  with  no  PML  is  used. 


IV.  Discussion 

A  PML  has  been  developed  for  applications  involving  ELFs 
in  lossy  media.  As  opposed  to  the  original  Berenger  PML  in  free 
space,  this  one  requires  a  split  JS-field,  but  not  a  split  H- field. 
Although  the  use  of  a  PML  is  not  as  crucial  as  it  might  be  in 


1081 

free  space  or  other  lossless  media,  it  substantially  decreases  the 
needed  computer  resources.  For  instance,  in  the  examples  in 
Section  IB,  it  was  found  that  a  problem  space  of  803  was  nec¬ 
essary  to  insure  that  the  boundaries  were  not  influencing  the  re¬ 
sults  when  there  was  no  PML.  After  the  PML  was  added,  the 
problem  space  was  reduced  to  60  x  20  x  20.  This  represents  a 
reduction  in  the  problems  space  from  512  000  cells  to  24  000 
cells. 
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Underwater  FDTD  Simulation  at  Extremely  Low 

Frequencies 

Yang  Xta  and  Dennis  M.  Sullivan,  Senior  Member,  IEEE 


Abstract — This  letter  describes  the  application  of  the  fi¬ 
nite-difference  time-domain  (FDTD)  method  to  the  simulation 
of  extremely  low  frequency  (ELF)  electromagnetic  signals  under 
water.  This  requires  substantial  modification  to  the  traditional 
FDTD  method,  as  well  as  the  development  of  an  analytic  method 
needed  to  verify  the  accuracy  of  the  FDTD  method. 

Index  Terms — Electromagnetic  propagation  in  absorbing  media, 
extremely  low  frequencies,  finite-difference  time-domain  (FDTD) 
methods. 


1.  Introduction 

MODERN  antiship  mines  can  be  detonated  by  the  electro¬ 
magnetic  signature  of  a  surface  ship  [  1  ].  For  this  reason, 
it  is  desirable  to  have  simulation  methods  to  study  the  propaga¬ 
tion  of  extremely  low-frequency  (ELF)  electromagnetic  waves 
under  water.  The  finite-difference  time-domain  (FDTD)  method 
[2],  [3]  is  one  of  the  most  widely  used  methods  in  electromag¬ 
netic  simulation.  However,  it  has  seen  limited  use  for  low  fre¬ 
quencies  in  lossy  media.  In  this  letter,  we  describe  the  use  of 
FDTD  for  ELF  simulation  under  water.  In  Section  II,  we  de¬ 
scribe  the  formulation  of  the  FDTD  method  that  has  been  found 
to  be  most  appropriate  for  this  application.  Section  III  describes 
the  method  of  two  equations,  two  unknowns  (2E2U)  that  is 
used  to  determine  the  resulting  amplitudes  when  the  FDTD  pro¬ 
gram  has  reached  steady  state.  Section  IV  describes  an  analytic 
method  that  was  developed  to  evaluate  the  accuracy  of  FDTD 
at  ELF.  Section  V  presents  an  example  of  ELF  simulation  in 
shallow  water,  the  type  of  problem  that  will  be  of  interest  for 
this  project.  Section  VI  ends  in  a  discussion,  including  remarks 
on  future  areas  of  research 


II.  Method 

We  begin  with  the  time-domain  Maxwell’s  equations 


dE 

r£°~dt  =  ^  X  H  ~ 

(la) 

on 

tL°~ot  =~VxE- 

(lb) 

Equation  (1)  is  for  three  dimensions,  but  for  the  purpose  of  il¬ 
lustration,  we  limit  the  discussion  to  the  Ex  and  Hy  fields  prop¬ 
agating  in  the  2  direction.  Equation  (la)  can  be  taken  into  the 
sampled  time  domain  using  the  usual  finite-differencing  proce¬ 
dures 

£?+!(*) -££(*)  (a\  +1 

"r  A t  +UJ  x 

1  //y  +1/2(fc  +  1/2)  -  //“ +l/2(fc  -  1/2) 
Co  At 

We  assume  the  cell  size  is  Ax  and  the  time  step  is  At.  The  £"  +  1 
can  now  be  calculated  from 


E?+1(k)  =  c a(k)E%(k) 

+c b(k)  [//yn+1/2(/c  +  1/2)  -  //;,+1/2(/c  -  1/2)]  (2a) 

where 


ca  = 


At  •  o 
€r£o 


- 1 


c  b  =  ca 


At 

(ere0Ax) 


(2b) 


There  is  a  crucial  choice  that  was  made  here.  Usually,  the  Ex 
term  next  to  the  conductivity  is  averaged  across  the  two  time 
steps 


£g+*(fc)-£g(fc)  +  fo_\  (E»+\k)  +  £"(*) 


A  t 


£o 


1  Hy+1/2(k  +  1/2)  -  IIy  +  1/2{k  -  1/2) 


Co  Ax 

which  would  lead  to  the  following  expression  for  ca: 


c  a  =  (  1  — 


At  •  a 
2>£r£  o 


1  + 


At  •  a 

2>£r£  o 


(3) 


At  ELF  frequencies  in  lossy  media,  the  ca  of  (3)  would  be  neg¬ 
ative,  leading  to  a  potentially  unstable  condition.  (The  imple¬ 
mentation  of  (lb)  into  FDTD  is  straight-forward  and  will  not  be 
presented  here). 

There  is  another  choice  that  leads  to  substantially  larger  time 
steps,  and  therefore,  substantially  faster  solutions  [4].  Once  the 
cell  size  Ax  is  chosen,  the  time  step  must  be  chosen  to  satisfy 
the  Courant  condition,  which  in  three  dimensions  is 
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At  < 


Ax 


-c 


max 


(4) 


where  cmax  is  usually  the  speed  of  light  in  a  vacuum.  The  mate¬ 
rials  that  will  be  of  interest  for  this  project  are  listed  in  Table  I 
(An  early  goal  of  this  project  is  to  study  propagation  in  lakes. 
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TABLE  I 

The  Properties  of  the  Materials  Used  in  the  Simulations  Described 
in  This  Paper  [5] 


Material 

£r 

cr  (SJm) 

Water 

80 

0.018 

Air 

1 

0 

Lake  water 

80 

0.018 

Mud 

40 

0.002 

Metal 

1 

107 

That  is  the  reason  that  lake  water  is  used  instead  of  sea  water.) 
The  complex  dielectric  constant  is  calculated  by 


cr  cr  '  • 

jve  o 

At  ELF  frequencies,  the  imaginary  part  of  the  dielectric  constant 
will  dominate  the  magnitude  for  all  the  materials  except  air. 
Therefore,  increasing  the  dielectric  constants  of  mud  or  metal 
to  80  would  make  very  little  difference.  If  we  assume  every 
material  in  Table  I  has  a  real  dielectric  constant  of  80,  then 
cmax  =  co/\/80  and  the  time  step  is  almost  an  order  of  magni¬ 
tude  greater,  (co  is  the  speed  of  light  in  a  vacuum).  Even  though 
air  is  one  of  the  materials  listed  in  Table  I,  air  is  a  boundary 
medium  in  this  project.  It  presents  almost  perfect  reflection  to 
an  electromagnetic  signal  in  water,  even  if  the  higher  dielectric 
constant  is  used. 

III.  The  Method  of  Two  Equations,  Two  Unkowns 

FDTD  is  a  time-domain  method.  Once  the  steady  state  has 
been  reached  for  a  simulation  problem,  it  is  desirable  to  know 
the  resulting  amplitude  and  phase  at  certain  locations  in  the 
problem  space.  For  frequencies  of  about  100  kHz  and  above, 
the  discrete  Fourier  transform  is  the  preferred  method  [6].  We 
have  found  that  at  ELF  frequencies,  the  method  of  two  equa¬ 
tions,  two  unknowns  (2E2U)  is  preferable  [7].  In  this  method, 
two  sample  points  are  taken 

E\  —  Asm(uJint\  -f  9)  (5a) 

E2  =  A  sin(u;in£2  4-  0).  (5b) 

Since  the  input  frequency  cvm  as  well  as  the  two  sample  points 
ti  and  t2  are  known,  the  only  unknowns  are  the  amplitude  A  and 
the  phase  9.  The  concept  of  solving  for  two  unknowns  from  the 
two  equations  is  straight-forward,  but  the  fact  that  the  inverse 
trigonometric  identities  must  be  taken  can  lead  to  inconsisten¬ 
cies.  It  has  been  found  expedient  to  add  an  offset  time 


to  each  of  the  times  1 1  and  t2 .  This  centers  the  two  sample  points 
symmetrically  on  the  ninety  degree  axis  and  avoids  problems 
when  taking  inverse  trigonometric  functions. 

IV.  Verification  of  the  Accuracy  of  the  Method 

An  analytic  solution  is  needed  to  verify  the  accuracy  of  the 
FDTD  method  at  ELF.  One  such  method  that  is  often  used  to 
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Fig.  I.  A  layered  dielectric  sphere  in  a  constant  E  field. 


verify  FDTD  formulations  is  a  layered  dielectric  sphere  illumi¬ 
nated  by  a  plane  wave.  A  Bessel  function  expansion  is  used  to 
calculate  the  resulting  fields  [8].  This  method  is  not  valid  below 
about  100  kHz. 

At  low  frequencies,  the  near  field  can  be  regarded  as  a  static 
field.  A  solution  for  a  layered  sphere  in  a  static  electric  field  was 
developed,  as  illustrated  in  Fig.  1.  A  spherical  boundary-value 
problem  has  solutions  of  the  form  [5] 

oo 

V(R,8)  =  £  [a„R"  4-  DUR-{1,+1)]  /J„(cos0).  (7) 

Pn(cos9)  are  the  Legendre  polynomials.  In  the  limit  far  from 
the  sphere,  V0(R,  #)|k->oo  —  —EoRcosO,  and  inside  the 
sphere  V\(R.9)  =  A\Rco^9.  There  are  two  boundary  con¬ 
ditions  at  a  dielectric  boundary:  E\t  =  E2t,eiEin  =  e2E2n. 
The  two  equations  resulting  from  the  boundary  conditions  are 

Bm  +  Amr®,  -  Dm+i  -  Am+ir-fn  =  0  (8a) 

and 

-  2-^Dm  +  —Amr?n  +  2 Bm+1 

£m  +  l  ^m  +  1 

-  ^m+ir,3,,  =  o.  (8b) 

The  constants  are  determined  by  Gaussian  elimination.  Once 
the  potential  V  is  known,  the  E  fields  are  determined  by 

which  can  be  converted  to  rectangular  coordinates. 

In  order  to  compare  the  FDTD  results  with  the  analytic 
method,  we  use  the  three-dimensional  problem  space  illus¬ 
trated  in  Fig.  2.  A  plane  wave  polarized  in  the  z  direction  is 
generated  at  one  end  and  subtracted  out  the  other  end.  The 
cells  used  in  the  simulations  are  five  meters  cubed  and  the  time 
steps  are  75  ns.  A  layered  sphere  with  dielectric  properties  to 
simulate  various  materials  lies  in  the  center  of  the  total  field. 
The  amplitude  of  the  Ez  field  is  determined  along  the  major 
axes  for  comparison  with  the  analytic  method  to  evaluate  the 
accuracy  of  the  FDTD  simulation.  These  axes  go  through  the 
sphere  and  extend  five  cells  out  in  each  direction.  The  problem 
space  is  50  cells  cubed. 

The  results  of  the  simulations  are  shown  in  Fig.  3  The  solid 
lines  are  the  analytic  results  and  the  circles  are  the  FDTD  values. 
Clearly,  the  comparisons  are  very  good.  In  Fig.  3(c),  there  is 
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Fig.  2.  The  configuralion  of  the  three-dimensional  simulation  space  used  to 
evaluate  the  accuracy  of  the  method. 


some  discrepancy  in  the  air  layer  in  the  middle.  This  is  not  too 
worrisome  because  in  our  problems  of  interest,  air  is  a  boundary, 
not  a  central  part  of  the  problem  space.  Note  that  this  condition 
was  not  caused  by  increasing  the  dielectric  to  80  as  discussed 
above;  the  same  result  is  obtained  using  a  dielectric  of  1  or  80. 


V.  Example 

Fig.  4  illustrates  the  type  of  simulation  of  interest  for  this 
project  Two  dipoles,  one  used  as  a  transmitter  and  one  as  a  re¬ 
ceiver,  are  submerged  in  shallow  water.  The  problem  space  is 
40  x  60  x  40  cells  and  each  cell  is  10  m  squared.  The  transmit¬ 
ting  current  is  simulated  by  the  H  fields  surrounding  the  middle 
of  the  transmitting  dipole 
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The  resulting  current  on  the  receiver  is  calculated  with  a  similar 
equation.  Each  simulation  required  2000  time  steps.  The  results 
are  shown  in  Fig.  5.  The  important  quantity,  i/,  the  transfer 
function,  is  the  ratio  of  received  current  to  transmitted  current, 
which  is  plotted  as  a  function  of  frequency. 


VI.  Discussion 

A  method  has  been  described  to  simulate  electromagnetic 
waves  propagating  under  water  at  extremely  low  frequencies. 
This  approach  necessitated  substantial  modification  to  the  usual 
FDTD  formulations.  Furthermore,  an  analytic  method  based  on 
the  Legendre  polynomials  was  developed  to  verify  the  accuracy 
of  the  FDTD  method. 

Those  familiar  with  FDTD  methods  will  notice  the  lack  of 
discussion  on  absorbing  boundary  conditions  (ABCs).  ABCs 
are  usually  required  to  prevent  outgoing  signals  from  being  re¬ 
flected  back  into  the  problem  space.  The  very  lossy  background 
medium  of  lake  water  has  prevented  this  from  being  a  concern 
for  the  examples  presented  in  this  letter.  However,  it  is  likely 
that  an  appropriate  ABC,  probably  one  based  on  the  perfectly 
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Xaxis  (m) 

(c) 

Fig.  3.  Comparison  of  the  FDTD  values  (circles)  versus  the  analytic  values 
(lines)  for  an  incident  plane  wave  at  I  kHz  and  a  layered  sphere  composed  of 
different  media,  (a)  The  inner  layer  is  mud,  the  outer  layer  is  mud  and  water 
(b)  The  inner  layer  is  metal,  the  outer  layer  is  mud.  (c)  The  inner  layer  is  air.  the 
outer  layer  is  mud. 


air 


Transmitter 


Receiver 


50  m 


200  m 

Fig,  4  Two  dipoles  are  submerged  in  shallow  water.  One  is  used  as  a  trans¬ 
mitter  and  the  other  as  a  receiver. 
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Frequency  (Hz) 

Fig.  5.  The  transfer  function  of  the  two-dipole  simulation  shown  in  Fig  4.  The 
frequency  range  is  3  Hz  to  3  kHz. 

matched  layer  (PML)  [9],  will  be  necessary  for  simulation  over 
long  distances. 

In  this  project,  it  is  anticipated  that  simulation  over  distances 
of  several  kilometers  will  be  required.  Some  form  of  near-to-far 
held  transformation  will  be  developed  to  model  the  EM  sources 
with  relatively  high  resolution  while  using  lower  resolution  to 
model  greater  distances  in  the  far  field 
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FDTD  Numerical  Tests  of  the  Convolutional — PML 
at  Extremely  Low  Frequencies 

Christopher  L.  Wagner  and  Jeffrey  L.  Young,  Fellow,  IEEE 


Abstract — Numerical  evaluation  of  the  finite-differencec 
time-domain  (FDTD)  convolutional  perfectly  matched  layer 
(CPML)  at  extremely  low  frequencies  (ELF)  is  conducted  herein 
to  arrive  at  acceptable  values  for  the  PML  parameters.  This  is 
accomplished  by  conducting  numerous  simulations  of  an  electric 
dipole  in  a  60  X  60  x  120  free-space  domain  and  by  benchmarking 
the  simulation  data  against  reference  data  for  strategic  observa¬ 
tion  points  within  the  domain.  Results  show  that  PML  attenuation 
on  the  order  of  60  to  70  dB  can  be  obtained  for  10  to  1000  Hz 
signals  in  the  quasi-static  region  of  the  dipole. 

Index  Terms — Absorbing  boundary  condition,  finite-difference 
time-domain  (FDTD)  methods,  perfectly  matched  layer  (PML). 


I.  Introduction 


IT  IS  well  known  in  oceanic  environments  that  only  ex¬ 
tremely  low-frequency  (ELF)  electromagnetic  waves  will 
propagate  over  long  distances  due  to  the  high  conductivity  of 
saltwater.  For  this  reason,  such  waves  are  quite  useful  in  com¬ 
munication  links,  or  can  he  undesirable  emissions,  as  caused  by 
high-powered  electric  drives  on  a  ship  platform.  In  either  case, 
the  propagation  characteristics  of  these  waves  can  be  understood 
from  computer  finite-difference  time-domain  (FDTD)  simula¬ 
tion,  particularly  in  the  littoral  region,  where  the  topological  and 
bathymetry  features  can  be  geometrically  complex.  To  accom¬ 
plish  such  a  simulation,  a  suitable  domain  truncation  technique 
is  needed  for  both  the  water  and  air  regions  of  the  domain.  Since 
the  ELF  signals  are  naturally  attenuated  in  the  water,  the  per¬ 
fectly  matched  layer  (PML)  development  for  the  air  is  the  most 
challenging. 

Classical  PMLs  used  in  FDTD  truncation  have  poor  per¬ 
formance  at  low  frequencies  and  potentially  suffer  late-time 
growth  [  1  ]— [3].  The  complex  frequency  stretching  scheme  in¬ 
troduced  by  Kuzuoglu  and  Mittra  [4]  alleviates  these  problems. 
The  FDTD  CPML  implementation  of  [4]  was  introduced  by 
Roden  and  Gedney  [5]  and  is  evaluated  here  for  ELF  perfor¬ 
mance.  For  this  work,  we  consider  1 0  to  1 000  Hz  to  be  the  ELF 
band. 
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II.  Formulation 

A.  The  FDTD  Problem  Statement 

The  FDTD  simulations  presented  herein  used  64 — bit  double 
precision  calculations.  All  simulations  are  performed  at  the 
Courant  stability  limit  to  minimize  dispersion  error  and  to 
advance  time  as  fast  as  possible.  The  FDTD  code  is  a  cubic  cell 
implementation,  with  a  cell  size  of  20  m.  The  medium  is  free 
space.  The  test  domain  size  is  60  x  60  x  120  cells,  including 
the  10-ccll-thick  PML.  The  electromagnetic  field  is  excited 
by  a  current  source  that  is  at  node  coordinate  (30,20,40).  A 
time -differentiated  Gaussian  waveform  is  used  as  the  excitation 
pulse.  This  pulse  has  no  dc  component,  so  no  persistent  charge 
will  be  deposited  into  the  grid,  which  would  produce  undesir¬ 
ably  large  dc  electric  fields  [6],  [7].  The  field  is  quantified  at  six 
observation  points  located  at  (30,10,80),  (30,20,80),  (30,30,80), 
(30,40,80),  (30,50,80),  and  (50,20,40),  respectively.  The  first 
set  of  grid  numbers  is  regarded  as  observation  point  I,  the 
second  set  as  observation  point  2,  etc.  Since  the  free-space 
wavelength  of  a  10-Hz  signal  is  3000  km,  it  is  clear  that  the 
observation  points  are  within  the  quasi-static  region  of  the 
dipole.  Such  near-field  observations  pose  significant  challenges 
to  FDTD  PML  development. 

Several  test  cases  are  considered  that  are  associated  with  var¬ 
ious  PML  parameters.  The  efficacy  of  each  PML  is  obtained  by 
benchmarking  the  FDTD  data  against  a  reference  solution,  as 
described  next. 

B.  Reference  Free-Space  Problem 

To  provide  a  reference  solution,  a  large  frec-space 
200  x  200  x  260  domain  with  perfect  electric  conductor 
(PEC)  walls  is  used.  The  geometry  of  the  source  and  receiver 
points  is  the  same  as  the  PML  test  cases,  but  the  free-space 
domain  is  larger  than  the  test  case  domains  by  140  cells  in  each 
direction.  This  reference  domain  is  large  enough  that  the  direct 
signal  is  fully  resolved  from  the  reflections  from  the  walls,  so 
the  reflections  can  be  removed  by  time-gating.  The  reference 
problem  only  needs  to  run  for  a  few  hundred  time-steps  to 
obtain  a  clean  direct  signal.  This  type  of  reference  solution 
includes  all  FDTD  numerical  errors,  thus  allowing  us  to  isolate 
the  effect  of  the  PML  induced  errors  from  all  others. 

C.  Signal  Processing 

To  extract  the  frequency  response  data  from  the  time-domain 
data,  fast  Fourier  transforms  (FFT)  are  used  [9].  The  simulations 
are  conducted  using  200  K  time  steps,  which  unfortunately  does 
not  give  sufficient  frequency  resolution  to  observe  the  ELF  re¬ 
sponse.  To  circumvent  this  problem,  the  time-domain  data  set 
is  extended  with  zeros  to  a  length  sufficient  to  obtain  the  lowest 
frequency  needed.  For  signals  that  decay  to  zero  (as  is  the  ease 
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with  good  PMLs),  zero  extension  is  proper.  For  numerical  pur¬ 
poses,  we  define  "zero"  as  less  then  10-12  relative  to  the  peak 
value.  The  zero  extended  data  sets  are  then  transformed  with  the 
FFT.  The  transformed  data  sets  are  then  used  to  compute  the  fre¬ 
quency-domain  performance  metric.  The  time-domain  and  fre¬ 
quency-domain  metrics  are  described  next. 

D.  The  Performance  Metrics 

To  measure  the  performance  of  the  PMLs,  an  energy  metric  is 
used.  The  energy  includes  all  the  field  components  in  the  metric. 
This  eliminates  the  possibility  of  choosing  an  especially  strong 
or  weak  field  component  at  random.  We  present  both  time-  and 
frequency-domain  metrics.  The  time-domain  metric  is  broad¬ 
band,  which  contains  all  spectral  information  contained  in  the 
excitation  signal.  The  frequency-domain  metrics  are  narrow- 
band,  calculated  at  selected  frequencies  of  interest. 

In  the  time  domain,  the  residual  energy  error  metric  is 

_  Yjt  {SweW  ±  6wh{t)) 

r_  E,  «(<)  +  <(<))  ’ 

where  w'e(t)  and  w'h(t)  are  the  electric  and  magnetic  energy 
densities  in  the  time-gated  reference  signal,  and  6uv(0>  6wh(t) 
are  the  residual  energy  densities  associated  with  the  PML.  The 
summations  are  over  the  full  simulation  time.  The  reference  en¬ 
ergy  electric  and  magnetic  densities  are  defined  as 

W'(0=^e'(<)-e'(<)  (2) 

and 

u-;,(o  =  i/ih  '(*)  •  h'(t)  o) 


where  e'(f)  is  the  time-gated  reference  FDTD  electric  field 
vector  and  h '(t)  is  the  reference  FDTD  magnetic  field  vector. 
The  residual  electric  energy  is  given  by 

bwe(t)  =  |e  (e(<)  -  e'(t))  •  (e(<)  -  c'(t))  (4) 

where  e  is  the  PML  FDTD  electric  field  vector  Similarly,  the 
residual  magnetic  energy  is 

Swh(t)  =  i/t  (h(t)  -  h'(t))  •  (MO  -  h'(0)  (5) 


where  again  the  primed  vector  is  the  reference  solution  and  the 
unprimed  vector  is  the  PML  FDTD  solution. 

In  the  frequency  domain,  the  residual  energy  error  metric  at 
angular  frequency  u>  is 


6W£(v)  +  6W„{u) 

W  eM  +  w„  M 


(6) 


where  the  residual  and  reference  energies  arc  defined  in  a 
manner  similar  to  the  time-domain  case. 


E.  The  PMLs 

In  the  frequency  domain,  the  CPML  tensor  coefficient  as 
given  by  Kuzuoglu  [4]  is 

£>‘j  —  +  Cf 1 1 —  X,y,Z.  (V) 


Frequency  (Hz) 


Fig.  I.  Hertzian  dipole  field  at  observation  poinl  1  for  several  different  PMLs 
with  /v  —  1 .  The  rapid  fall-off  of  amplitude  afler  I  MHz  is  due  to  the  limited 
bandwidth  of  lhc  source  current  waveform. 


The  real  coordinate  stretch  k,  conductivity  <x,  complex  fre¬ 
quency  stretch  a,  and  their  polynomial  scaling  characterize  the 
PML. 

The  test  PMLs  are  10  cells  thick,  with  the  parameters  having 
polynomial  scaling.  The  conductivity  a  and  coordinate  stretch 
k  use  a  fourth-order  scaling  polynomial,  while  the  complex  fre¬ 
quency  stretch  a  uses  third-order.  As  is  usual  with  the  CPML, 
the  scaling  polynomial  for  the  conductivity  and  real  coordinate 
stretch  increases  into  the  PML,  while  the  complex  frequency 
stretch  polynomial  decreases  into  to  PML.  The  maximum  con 
ductivity  is  set  according  to  the  optimum  [8]  given  by 

771  +  1 

where  m  is  the  polynomial  order,  Ax  is  the  space  grid  size,  and 
er  is  the  relative  dielectric  constant  The  maximum  complex 
frequency  stretch  a  is  set  by 

=  2jne0fa  (9) 

where  fa  is  the  CPML  break  frequency. 

The  problem  is  to  hnd  ranges  for  the  parameters  fa  and  k  that 
provide  good  performance  at  ELF.  This  can  only  be  done  em¬ 
pirically.  Representative  test  cases  are  provided  next  to  demon 
strate  this  empirical  process. 

III.  Results 

To  validate  the  simulations,  the  exact  Hertzian  dipole,  fre¬ 
quency-domain  solution  is  compared  to  the  transformed  FDTD 
simulation  data  in  Fig.  1 .  The  field  is  observed  at  point  1 .  When 
fa  is  small,  the  low-frequency  performance  of  the  PML  is  poor. 
Likewise,  when  fa  is  excessively  large,  the  PML  performs 
poorly  at  high  frequencies.  A  PML  with  a  reasonable  value  of 
the  break  frequency  provides  a  simulation  that  closely  matches 
the  theoretical  prediction  over  the  full  excited  frequency  band. 
We  have  found  empirically  that  fa  «  1  MHz  seems  optimum 
for  ELF  simulations.  This  conclusion  is  also  valid  when  the 
field  is  observed  at  other  strategic  observation  points,  i.c.  2,  3, 
4,  5,  and  6. 
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Fig.  2.  Time-domain  results  at  observation  point  1  for  k  =  1.  For  limes  less 
than  about  8c-6  s,  the  direct  signal  is  seen.  After  the  direct  pulse  has  passed  by 
the  sample  point,  various  levels  of  residual  fields  are  seen. 


Fig.  3.  Normalized  PML  frequency  response  at  observation  point  1  for  k  =  1 . 
Lower  amplitudes  are  better  performing  PMLs. 

See  Fig.  2  for  the  time-domain  performance  of  four  test  and 
the  reference  simulations.  The  reference  simulation  shows  the 
direct  signal  clearly  separated  from  the  reflected  signals,  the 
latter  of  which  can  be  removed  by  time-gating,  as  noted  previ¬ 
ously.  With  fa  —  1  Hz,  there  is  slow  long-term  decay.  (In  some 
simulations,  when  fa  —  0,  slow  growth  has  even  been  reported 
[3].)  As  fa  is  increased  to  the  optimum,  the  absorption  increases 
relative  to  the  fa  ~  1  Hz  case.  As  fa  is  increased  further,  the 
absorption  degrades,  but  still  with  good  late-time  fall-off. 

In  Fig.  3,  the  frequency-domain  residual  energy  error  metric 
as  computed  by  (6)  is  shown  for  observation  point  I .  In  this  plot, 
a  better  PML  will  have  a  lower  response.  The  PML  is  tested 
with  various  break  frequencies  /a,  each  with  n  =  1  With  fa 
too  large  or  small,  there  is  poor  PML  absorption.  As  can  be  seen 
in  both  the  time-domain  and  frequency-domain  plots,  there  is  an 
optimum  value  for  the  break  frequency  for  ELF  simulations.  For 
the  tests  performed  here,  fa  ^  1  MHz  provides  the  best  PML 
absorption  at  ELF  with  a  relative  error  on  the  order  of  0.02%. 

There  is  up  to  a  factor  of  100  variation  in  error  in  PMLabsorp 
tion  across  the  six  sample  points,  as  shown  in  Fig.  4.  Surpris¬ 
ingly,  sample  point  I  has  better  ELF  performance  than  sample 
point  5.  From  a  wave  perspective,  sample  point  1  is  the  grazing 
incidence  case;  however,  given  that  the  fields  are  quasi-static. 


Fig.  4  PML  frequency  response  for  fa  =  10°  Hz  and  *  =  1  for  the  six 
observation  points. 


K 


Fig.  5.  Time-domain  PML  residual  energy  error  metric  at  observation  point  1 
for  various  values  of  the  real  coordinate  stretch  k  and  the  break  frequency  fa . 
FDTD  simulations  arc  performed  at  each  gnd-line  intersection 

grazing  incidence  has  no  real  meaning.  Clearly  from  Figs.  3  and 
4,  good  choices  of  parameters  provide  better  than  70  dB  of  PML 
absorption  over  the  ELF  band  at  favorable  observation  locations 
and  more  than  60  dB  attenuation  in  unfavorable  locations. 

A.  Time-Domain  ( Broadband )  Performance 

Figs.  5  and  6  show  the  time-domain  contours  of  10logi0r 
as  k  varies  from  1  to  20  and  as  fa  varies  from  1  to  I010  Hz 
for  observation  points  I  and  5,  respectively.  There  is  variation 
in  the  location  and  depth  of  the  global  minimum  across  the  six 
observation  points.  However,  with  fa  «  1  MHz  and  3  <  k  < 
8,  the  PML  provides  —70  dB  or  better  performance  at  all  six 
observation  locations. 

B.  Frequency- Domain  ( Narrowband )  Performance 

Figs.  7-9  show  the  frequency-domain  PML  performance  at 
100  Hz,  10  kHz,  and  1  MHz,  respectively.  The  1  Hz  to  I  kHz 
optima  vary  widely  with  k,  /a,  and  observation  position  For 
some  specific  lest  locations,  frequencies,  and  PML  parameters, 
there  are  very  deep  minima  on  the  order  of  —  120  dB  in  some 
cases.  Because  the  location  and  parameters  of  these  minima  do 
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Fig.  6.  Time-domain  PML  residual  energy  error  metric  at  observation  point  5. 


* 


Fig.  7.  100-Hz  frequency-domain  PML  residual  energy  error  metric  at  obser¬ 
vation  point  I.  At  the  same  observation  point,  the  1  to  100  Hz  results  look  es¬ 
sentially  identical. 

not  vary  in  a  regular  way,  it  is  important  to  select  operating  con¬ 
ditions  for  the  PML  based  on  the  performance  at  several  obser¬ 
vation  locations.  For  example,  at  10  Hz,  100  Hz,  and  1  kHz, 
fa  zz  10°  Hz  with  1  <  k  <  20  provides  —60  dB  or  better  PML 
performance  at  all  six  observation  points.  If  18  <  n  <  20,  then 
the  PML  provides  —70  dB  or  better  performance  at  all  six  ob¬ 
servation  points.  Apparently,  the  narrowband  metrics  at  ELF  are 
improved  with  larger  values  of  k  as  compared  to  the  wideband 
time-domain  metric.  Fig.  9  shows  that  at  1  MHz,  smaller  values 
for  k  provide  the  best  PML  absorption.  The  optimum  PML  pa¬ 
rameters  depend  on  the  metric  used  and  on  the  frequencies  of 
interest, 

IV.  Conclusion 

For  ELF  PML  development,  our  empirical  research  shows 
that  when  fa  zz  I  MHz,  n  ^  6,  a  —  aopt,  fourth-order  polyno¬ 
mials  for  a  and  n  are  invoked,  and  a  third-order  polynomial  for 
a  is  invoked,  then  at  least  60  dB  of  PML  attenuation  is  obtained 
in  the  quasi-static  region  (i.e.  very  near-field)  of  the  dipole,  for 
both  the  wideband  and  ELF  narrowband  metrics.  This  is  more 
than  adequate  for  high-quality  FDTD  simulations  in  the  ELF 
band. 
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Fig  8  10-kHz  frequency-domain  PML  residual  energy  error  metric  at  obser¬ 
vation  point  1 . 
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Fig.  9.  I -MHz  frequency-domain  PML  residual  energy  error  metric  at  obser¬ 
vation  point  I 
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Abstract — The  propagation  of  electromagnetic  fields  from  an 
electric  or  magnetic  dipole  near  an  interface  has  been  described 
mathematically  for  man)  years  using  Sommcrfeld  integrals.  In 
the  high  frequency  case  for  which  the  dipole  and  field  point  are 
buried  in  the  lower  (and  much  lossier)  medium,  simple 
approximations  have  been  derived  that  can  be  interpreted  as 
up-over-and-down  propagation  (i.e.,  plane  wave  from  source  to 
the  surface  -  surface  wave  to  a  point  above  the  field  point  -  plane 
wave  to  the  field  point).  In  this  paper,  a  simple  expression  \alid 
for  the  low  frequency  case  (i.e.,  relevant  dimensions  in  the  upper 
medium  are  electrically  small)  is  derived.  In  this  case  the 
up-over-and-down  propagation  is  described  as  near  field 
propagation  from  source  to  the  surface,  quasi-static  propagation 
along  the  interface  to  above  the  Held  point  and  plane  wave 
propagation  to  the  field  point. 

Index  Terms — Electromagnetic  fields,  electromagnetic 

propagation,  extremely  low  frequency. 

I.  INTRODUCTION 


good  approximations  to  the  electric  and  magnetic  fields  are 
obtained  Based  on  these  approximations,  the 
up-over-and-down  behavior  is  observed  and  discussed.  The 
dipole  source  here  is  chosen  to  be  a  horizontal  electric  dipole 
(HED).  The  HED  was  selected  because,  using  achievable 
dipole  moments  and  commonly  available  receiving  equipment, 
it  can  be  shown  that  the  HED  fields  are  detectable  at  larger 
distances  than  those  of  other  dipole  types  (i.e  ,  vertical  electric 
dipole  (VED)  or  vertical  (VMD)  or  horizontal  (HMD) 
magnetic  dipoles).  To  be  more  specified,  it  was  assumed  from 
this  study  that  the  maximum  dipole  moments  for  electric  and 
magnetic  dipoles  are  50  A-m  and  2500  A-m  respectively  and 
that  the  minimum  detectable  electric  and  magnetic  fields  are 
lgV/m  and  40pA/m  respectively.  Using  these  values,  the 
horizontal  electric  field  component  that  is  perpendicular  to  the 
HED  direction  can  be  detected  to  a  distance  of  800  meters  to  the 
source.  No  other  field  component  from  any  other  dipole  can  be 
detected  beyond  about  200  meters. 


THE  electric  (E)  and  magnetic  (H)  fields  due  to  a  dipole 
(vertical  or  horizontal,  electric  or  magnetic)  buried  in  a 
conducting  half-space  have  been  well  studied  for  decades 
and  can  be  written  in  terms  of  Sommerfeld  integrals  [1]  -  [4] 
When  the  dipole  and  the  observation  point  are  both  in  the 
conducting  medium  close  to  the  interface  relative  to  their 
horizontal  spacing  and  the  frequency  is  “low”,  it  is  possible  to 
interpret  the  propagation  mechanism  as  a  simple 
up-over-and-down  process.  Here,  up-over-and-down  means 
that  the  field  propagates  vertically  up  crossing  the  interface  to 
the  free  space  medium,  then  propagates  horizontally  along  the 
interface,  and  finally  propagates  vertically  down  to  the 
observation  point.  While  this  behavior  is  somewhat  similar  to 
the  high  frequency  phenomenon  observed  by  previous  authors, 
it  is  also  different  because  the  fields  in  the  free  space  region  are 
quasi-static  [5]. 

The  formulas  for  the  fields  given  by  Sommerfeld  integrals 
are,  while  exact,  very  complicated.  In  this  paper,  the  integrals 
are  simplified  by  using  some  reasonable  assumptions  given  the 
range  of  parameters  of  interest.  Then  a  set  of  simple  but  very 
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II.  Geometry 


The  geometry  of  the  model  is  shown  in  Fig.  LA'/  oriented 
HED,  which  has  a  dipole  moment  of  Idl,  is  on  the  V  axis  and 
buried  ‘/T  meters  below  the  surface  in  a  conducting  half-space 
(z  <  0).  The  upper  half  space  (i.e.,  z  >  0)  is  assumed  to  be  free 
space. 


Free  space  (#0) 


Z 

z  =  0 


co-  cr0,  /f0 


observation 

HED4  l  =  -h  P°'nt 


Conducting  medium  (#1)  Cv  CT,,  fjQ 


Fig  1  Geometry  of  the  model 

As  noted  in  the  figure,  c,  and  a,  are  the  permittivity  and 
conductivity  of  the  ih  half  space  (/  =  0  and  1  for  free  space  and 
conductor,  respectively),  e,  -  £n£0,  where  £n  is  the  relative 
permittivity  and  £0  is  the  permittivity  of  free  space.  It  is 
assumed  that  all  materials  have  the  permeability  of  free  space 
/i0.  The  cylindrical  coordinate  system  (p,  (p,  z)  is  used  in  this 
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paper,  where  x  -  pcostp  and  y  =  psin^. 


111.  SOMMFRFELD  INTEGRAL  METHOD 

The  Sommerfeld  integral  method  will  be  very  briefly 
introduced  here.  This  method  uses  the  integral  representations 
of  vector  potentials  to  determine  the  E  and  H  fields.  To  find  the 
E  and  H  fields  due  to  the  HED,  two  non-zero  components  of 
vector  potential  are  required  [6,  7].  Here,  the  y  and  z 
components  of  the  magnetic  vector  potential,  Ay  and  AZi  are 
chosen. 

K  =  K,  [fx{X)e-^XJ0{Xp)dX  (z  >0)  (1) 
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where  JQ(Xp)  is  the  Bessel  function  of  the  first  kind  of  order  zero 
and 
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s\  -  sf  —  j  cjco  is  the  complex  permittivity  of  the  i  half  space 

(/  =  0  or  1 ),  kt  is  the  wave  number  where  Re(^)^0  and  Re(w,)  ^0 
defines  the  proper  Reimann  sheet  of  the  complex  plane.  The 
first  term  in  (2)  is  the  source  term  and  R  -(p2  +  z2)]f2  is  the 
distance  from  the  dipole  to  the  observation  point.  lyJ  and  l:i 
represent  the  integral  terms  in  (2)  and  (4),  respectively.  The 


source  term  K 


R 


in  A1  is  the  vector  potential  of  the 


dipole  itself  in  an  infinite  homogeneous  conducting  medium  It 
can  be  written  in  integral  form  as 
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Functions  g\  and  g2  are  arbitrary  coefficient  functions  of 
the  integration  variable  2.  They  are  determined  by  matching  the 
boundary  conditions  at  z  =  0.  Given  the  vector  potentials  Ai 

and  A ,  the  E  and  H  field  can  be  obtained  from: 
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By  matching  boundary  conditions,  i.c.,  the  tangential 
components  of  E  and  H  fields  are  continuous  cross  the  z  =  0 
plane,  the  coefficient  functions  can  be  determined  as 
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The  E  and  H  fields  in  the  conducting  medium  can  be 
formulated  by  inserting  (2),  (4),  ( 1 3)  and  ( 1 4)  into  (6)  -  ( 1 1 ). 


IV.  Simplification  of  the  integral  for  a\ 

The  objective  of  this  paper  is  to  derive  simple  but  acceptable 
approximations  for  the  fields,  which  can  be  interpreted  to 
provide  good  insight  into  the  physical  behavior  of  the  wave 
propagating  from  source  to  receiver.  One  fundamental  problem 
with  evaluating  the  integrals  shown  in  ( 1 )  to  (4)  is  that,  for  large 
values  of  p  compared  to  h  and  z,  the  rapid  oscillations  of  the 
Bessel  function  cause  difficulties  with  the  numerical 
integration.  To  remedy  this  problem,  the  contours  of  integration 
will  be  deformed  in  the  complex  plane  so  that  the  integrand 
decays  exponentially  for  large  values  of/).  This  transformation 
will  also  allow  other  simplifying  approximations  that  will  lead 
to  a  simple  interpretation  of  the  final  result. 


A  Deformation  of  The  Integral  Contour 

If  (14)  is  inserted  into  the  integral  portion  of  (4)  and  the 
exponential  term  is  removed  from  g2(X)y  this  integral  becomes 
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(16) 


Using  the  identities 
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ioW^[^"W+w'21w] 

H"\-x)  =  -H(02\x) 

where  Hq\x)  and  H{02)(x)  are  the  Hankel  functions  of  the 

first  and  second  kind  of  order  zero,  respectively,  the  integral 
range  in  (15)  can  be  expanded  to  (-00,  +00).  Since  w0,  u  1,  and 
g'  (A)  are  all  even  functions  of  X 

4c«;we“i<!‘X,(^  (|7) 

For  the  function  g'2(X)>  there  is  one  pole,  Xp,  and  two  branch 

points,  ko  and  k\9  in  the  complex  X  plane.  The  branch  cuts  are 
selected  to  be  vertical  lines  from  the  branch  points  to  negative 
infinity.  Then  the  integral  contour  in  (17)  can  be  deformed  into 
a  contour  CB  which  is  illustrated  with  the  dashed  line  in  Fig.  2. 
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Fig.  2.  Deformation  of  the  integral  contour 

With  this  deformation  and  the  fact  that  the  Hankel  function 
goes  to  zero  exponentially  along  the  infinite  semi-circle,  the 
integral  along  the  real  axis  is  converted  to  the  residue  of  the 
pole,  R^,  plus  the  integrations  along  C\  — >  C4,  which 
encompass  the  two  branch  cuts.  Thus 

Z  *m+c  a+t-j +( * 

For  | A:, |  »  |^o|,  the  integral  along  the  branch  cut  of  k\  is  much 
smaller  than  that  of  k$  and  can  be  ignored.  Note  that,  in  the 
complex  plane,  the  sign  of  u0  will  change  when  crossing  the 
branch  cut  associated  with  k 0.  Given  the  choice  of  branch  cut, 
Re(wo)<0  and  Re(wo)>0  on  the  left  and  right  sides,  respectively, 
as  shown  in  Fig.  2.  In  addition  while  the  pole  is  in  the  proximity 
of  the  branch  cut  integration  and  is  evident  in  the  integrand,  its 
contribution  to  the  integral  is  negligible  for  the  low  frequencies 
considered  here.  Thus  the  pole  residue  can  be  ignored. 
Therefore 


/„=4f.  g'M)e^-h)XH™(Xp)JX  (18) 

It  has  been  shown  that  ( 1 8)  is  valid  when  /?,  z  «  p ,  l*.l  »  |Aol 
and  \k\p\  »  1 . 

The  parameters  and  their  values  /  value  ranges  used  for  the 
simulations  are  listed  in  Table  1.  The  conductivity  and 
permittivity  of  the  lower  medium  represent  typical  lake  water. 
They  will  also  be  used  for  all  the  following  simulations  in  this 
paper. 


TABLE  1 


List  of  parameters  and  their  values  /  value  ranges 


Free  space 

Conducling  medium" 

Relative  permittivity, 

(e,  =  0) 

1 

1 

Conductivity,  <7  (S/m) 

0 

0018 

Permeability,  n  (H/m) 

4it*10‘7 

4**10-’ 

Dipole  depth,  h(m) 

20 

z(m) 

-10 

Horizontal  distance,/? 

i"»i  . 

100-  10000 

Dipole  moment,  tdl 
(A-m) 

1 

Frcquency,/(Hz) 

100-  3000 

B.  Simplification  of  The  Integrand 


Since  for  p  □  |  /;  |,  |  z  | ,  the  decay  of  the  integrand  along  C| 
and  C2  is  controlled  by  the  value  of  \Xp\,  the  integral  can  be 
truncated  at  |Ap|=10  and  since  we  assume  |2p|»l 

Ux=(X2-k?)U2*jkx  (19) 

With  (19)  the  exponential  term  in  (18)  can  be  extracted  from 

the  integral,  which  leads  to 

Jk 

4  *  — —  f.  g2{X)xn{2\xP)dx  (20) 

Z  ^  l+'  2 

Further  since  (i.e.,  quasi-static  forz>0),  it  is  reasonable 

to  assume  that 

u0  =  ±(X2-k02),/2  *±X  (21) 


because  \ko\  is  very  small  compared  to  \X\  over  the  largest 
portion  of  the  integral.  The  approximations  have  been  made 
here  can  be  summarized  as 

pti  | H  \z\ 

\koP\]  1 

l*,pf  1 


Now,  if  the  approximations  (19),  (21)  and  |e'i|»|£o|  are 
made  in  (16)  then 
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lf  8?t)W  =  "32  2-,  7  and  g'P(X)  =  —  2  —  ,  represent 
A  +  JKXA  A  -jk^A 

g'2(A)  on  the  right  side  and  left  side  of  the  branch  cut  of  £0, 

respectively,  then  the  integral  in  (20)  can  be  approximated  as 

Ki  “  l  7 A—H<2\Xp)dX+l  -l.-Hi2\Xp)dX 
*1  A-jkx  *  2  A  +  jk , 


(22) 


Using  the  asymptotic  approximation 

H{02)(Ap)zz  for  large  \Ap\  (i.e.,  most  of  the 

\  nAp 


integral) 
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(23) 


4l{X+jkx) 

Letting  X  =  k$  -  js ,  and  changing  the  integral  variable  from  X  to  s , 
/•i  finally  becomes 

[27  [  r  je~p'ds 
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Ignoring  &o  in  both  the  denominators  of  the  two  integrands, 
which  is  reasonable  because  \ko\  is  very  small  compared  to  \s\ 
over  most  of  the  integral,  /  2]  reduces  to 


where 


'.-r 


'.-r 


yfs  (s-k^ 


(25) 
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Using  this  result  in  /ls  results  in 


;  Vp) 


Similarly  for/2 


a  ~-'Tx/{k\4p)- 

Therefore,  the  integral  in  (25)  is 

/;  *_4V2y— — 

Ap 

The  integral  /rl  then  can  be  approximated  as 
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(26) 


It  can  be  shown  that  (26)  is  approximately  a  factor  of  1.4 
larger  than  the  exact  result  in  (15)  and  that  this  difference  is 
relatively  stable  over  a  wide  range  of  parameters.  Given  this 


and  the  fact  that  an  attempt  to  find  a  missing  >/2  factor  did  not 
succeed,  a  further  study  of  the  approximation  used  to  derive 
(26)  was  carried  out.  This  study  indicated  that  the  dominant 
part  of  the  error  resulted  from  the  replacement  of  the  Hankel 
function  by  its  asymptotic  expansion.  Given  this,  a  correction 
term  can  be  written  as 
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ds 


Clearly  most  of  the  contribution  to  this  integral  comes  from 
small  values  of ps.  Thus  the  integral  is  (somewhat  arbitrarily) 
truncated  at  ps  -  B  =  0.3  and  the  Hankel  function  is  replaced  by 
its  small  argument  expansion.  Given  this,  the  correction  term  is 
written  as 
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ds 


Therefore,  the  calculation  of  the  complicated  integral  in  (18)  is 
reduced  to  the  problem  of  evaluating  the  two  relatively  simple 
integrals  in  (25).  I\  can  be  analytically  evaluated  as  [8] 

where  D_x  ^2kxp  j  is  the  parabolic  cylinder  function  of 
argument  ^2 kxp  with 

The  second  term  in  the  bracket is  the  probability 
integral,  which  has  the  asymptotic  approximation  for  |£,/>|»I: 


where  inside  the  integral  e,k"p  □  e~p*  □  1  and  M-B/p.  The 
integral  in  (28)  can  be  analytically  evaluated  and  the  result  is 

(29) 
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let  and  /('3  can  be  further  simplified  by  expanding  the  natural 
logarithm  function  in  Taylor  series. 

T"  jM 


4-1  ^-2c~e-^  and  /„  *  8J- 
k,P  V 
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The  correction  term  is  then  rewritten  as 
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(30) 

is  a  constant  It  is 


interesting  to  note  that  the  functional  dependence  of  (30)  is 
identical  to  that  of  (26).  Hence  adding  (30)  to  (26)  results  in 
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/,i« 


eMz-h)  je~JkoP 
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(3D 

2  kxp 

which  is  identical  to  (26)  except  for  the  constant  and  that  this 
constant  is  approximately  1/1 A  times  the  constant  in  (26)  when 
M-  0.3  (i.e.,  =  1 .72).  Then  (31)  becomes 

ry-J^P 


/zl  ~  -2  jejk'u~>,)  ■  ■ 


k,P 


(32) 


When  /?,  z  «  p,  | k\\  » |/^|  and  \kxp\  »  1 ,  (32)  approximates  the 
exact  integral  very  well.  These  conditions  are  roughly  mapped 
to  the  following  range  of  parameters:  hyz<  100m,  100Hz  < /< 
3000Hz,  500m  <  p  <  1 0000m  and  0.00 1  S/m  <  a  <  1  OOS/m.  The 
error  of  (32)  compared  to  the  exact  integral  of  lz\  in  (15)  is  less 
than  10%  when  100Hz  < /<  3000Hz  and  500m  <  p  <  10000m. 
Fig.  3  shows  the  comparisons  of  magnitude  and  phase  angle 
between  the  approximation  in  (32)  and  the  exact  integral  in 
(15),  with  dipole  frequency  of  1000Hz. 
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results  for  the  E  and  H  fields.  Therefore,  (32)  will  be  used  as  the 
approximation  of  IzX. 

V.  Simplification  of  the  integral  for  a\ 

If  the  exponential  term  is  pulled  out  from  f2(X)t  Iyl  of  (2)  can 
be  rewritten  as 

lyX=[nW'U-k)MMp)dl  (33) 

where  /2'(A)=- 


u\  (uo  +  u\) 

At  this  point,  the  strategy  for  simplifying  ly\  is  the  same  as 
that  used  for  IzX.  First,  the  integral  contour  is  expanded  and 
deformed.  The  deformed  contour  is  the  same  as  that  illustrated 
in  Fig.  2  except  that  there  is  no  pole  in  this  case. 

/„=M  P(X)e"*-h'XH™{Xp)dX 

Z  *  B 

Again,  using  the  argument  that  the  integral  along  the  branch  cut 
of  ko  dominates  the  total  integral,  IyX  can  be  approximated  as  the 
sum  of  the  integrals  along  C\  and  C2 

,  PW^-h)XH^\Xp)dX  (34) 

2  ^  1 i 

The  approximations  given  in  (19)  and  (21)  still  work  and 

Ik  2 

given  these  the  function  reduces  to  /2f(A)  ~ — — - 

A{A±x) 

The  change  of  signs  in  it  is  due  to  that  w0  takes  different  signs 
on  the  left  and  right  sides  of  the  branch  cut  of  k Q.  The  integral  in 
(34)  becomes 


Iyl*  2  A 


f  XH?>(Xp)dX 

*■  1  ik  —  A 


A  ~  x 


+  f  — -7  XH™ (Xp)d X 

*  2  jkx  +  A 

Combining  the  two  integrals  on  the  right  hand  side  results  in 

r  4  ;2 

Use  the  asymptotic  approximation  of  the  Hankel  function 

H{n2)  ( Ap)  *  1-^-  ■  e~)Xfi ,  and  the  variable  change  X  =  ko  *  js, 
0  \nkp 

ly]  becomes 

/  .  «  -2 j  fZLe,l'u-h>e-jt°p  ■  r-&zM!Le-~d, 

\trp  ■*)  k[  +(k0-js ) 

Since  |5|»|  ko  |,  all  ko'  sin  the  integrand  can  be  ignored  and 
~2 


Fig  3  Exact  integral  of  lti  vs  approximation,  f-  1000Hz  (a)  magnitude,  (b) 
phase  angle 

This  analysis  is  helpful  for  understanding  the  error  incurred 
during  the  derivation  of  the  approximation  for  (15).  Further,  the 
correction  term  significantly  reduces  the  error  and  can  be  easily 
calculated.  It  is  shown  in  (32)  and  can  be  used  to  derive  simple 


‘yi 

3/2 


I  np 


(35) 


where  /,  =  f*  — - r-e  p'ds  .  From  the  table  of  integrals 

3  *52-*,2 

[8],  /3  can  be  analytically  evaluated  and  further  simplified  as 
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3yfn  1 


4^i  P 


5/2 


and  the  integral  lyl  is 


1  ~  e  >k"P 


2k 


(36) 


P 

Numerical  calculations  indicate  that  the  approximation,  (36), 
is  approximately  5%  larger  in  magnitude  than  the  exact  integral, 
(33).  This  error  is  relatively  stable  over  the  parameter  range 
100Hz  <  f  <  3000Hz  and  p  >  500m.  Since  the  error  is  small  for 
this  case,  it  is  not  necessary  to  add  a  correction  term  to  (36). 
Rather  the  factor  of  3>/2/4al  06  is  simply  set  equal  to  1 
resulting  in: 


/  _ _ 

~  k; 


~JkoP 


(37) 


Again,  when  h ,  z  «  py  |A||  »  \kj(  and  |Aip|  »  1,  (37) 
approximates  the  exact  integral  (33)  very  well.  Fig.  4  (a)  and 
(b)  give  the  comparisons  of  the  magnitude  and  the  phase  angle, 
respectively,  between  the  exact  integral  of  (33)  and  its 
approximation  (37).  For  this  case,  the  magnitude  error  of  (37)  is 
even  less  than  6%  when  p  >  500m. 
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Fig  4  Exact  integral  of  Iy,  vs  approximation  (a)  magnitude,  (b)  phase  angle 


VI.  The  approximations  for  E  &  H  field 

Given  the  approximations  (32)  and  (37)  for  the  integral 
portions  of  the  vector  potentials  A\  and/^,  respectively,  the 

electromagnetic  Fields  can  be  found.  The  approximate  results 
for  A |  and  Iyl  in  (32)  and  (37),  respectively,  can  be  used  in  (2) 
and  (4)  and  (6)  ~(1 1)  to  find  the  complete  set  of  E  and  H  fields 
in  region  #1.  The  components  of  total  E  and  H  fields  can  be 
written  as  the  combination  of  the  source  terms  and  the  reflected 
terms. 

E\  =  El  +  E[-  E'y  =  El  +  E'>r-  E]  =  El +  E't/, 

K  =  Hi  +  Hi;  H\  =  Hi  +  H\r;  H\  =  H'm  +  H\r. 
where  the  components  with  * s'  in  the  subscript  refer  to  the 
source  terms  of  the  fields  and  those  with  V  in  the  subscript 
refer  to  the  reflected  fields.  The  source  terms  of  the  fields  arc 

-Jk\R 

=  Ax'~ — -3yA:|/?-3j-sin^cos^  (38a) 

R 

El=Ae-^-[(-ktR2  +  jk,R+\) 

R  LV  7  (38b) 

+  (£2/?2-3;A:!/?-3)(sin^)2  ] 

El  =  A  •  +  7  {k;R2  - 3  jk,R- 3)  sin  <f>  (38c) 


Hi  = 


Idl  (z  +  h)e~,k'R 


At: 


Ri 


(jk,R  + 1)  (38d) 

(38c) 

Jr] I  P~jk'R 

H'a  -r-{jkxR  +  \)cos<f>  (380 

4  71  a 


//;=o 


where  A,  = 


jldl 

^no)e\ 


.  The  simplified  reflected  fields  are 


El*2jAXeJ*'u-h)A 
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F'^-ljAX-e 


„-)KsP 
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The  results  for  E  and  H  fields  calculated  by  (39)  are 
compared  to  the  results  from  the  exact  Sommerfeld  integral 
formulas.  Fig.  5  shows  the  comparing  result  of  the  reflected 
field  El  •  For  the  reflected  magnetic  field,  similar  result  can  be 

obtained. 


P(m) 

Fig  5.  Exact  reflected  field  Eyn  vs  its  approximation  (39b),  0  =  k 


VII.  Discussion 

At  any  observation  point  P,  (p,  <p,  z),  in  the  lower  conducting 
medium  (z<0),  the  total  field  is  the  combination  of  the  incident 
field  directly  from  the  dipole  and  the  reflected  field  due  to  the 
interface,  as  shown  in  (38)  and  (39),  respectively.  In  the  far 
field  of  the  conducting  medium  (|A|p|»/?,  z),  the  incident  field 
decays  much  faster  than  the  reflected  field,  which  suggests  that 
the  propagation  mechanism  may  involve  some  fields 
propagating  in  the  low  loss  free  space  region.  For  the  set  of 
parameters  considered  the  source  term  can  be  ignored  beyond 
approximately  500m. 

In  fact,  the  approximations  (39a)  ~  (39f)  can  be  interpreted  to 
have  an  up-over-and-down  behavior,  similar  (but  not  identical) 
to  the  mechanism  that  has  been  studied  in  propagation  of  high 
frequency  radio  waves  near  a  boundary  [5].  The  extremely  low 
frequency  (ELF)  case  is  different  because  the  fields  in  air  are 
quasi-static.  Here,  some  insight  into  this  behavior  will  be  given. 

Consider  the  Ex  component,  which  is  perpendicular  to  the 
HED’s  orientation  and  easier  to  be  detected,  at  the  observation 
point  P  (p.  <p,  z).  As  already  discussed,  the  total  field  component 
can  be  approximated  by  the  corresponding  reflected  field 
component  in  the  far  field.  Thus,  the  total  x  component  of  the 
electric  field  can  be  approximated  by  (39a).  In  the  far  field  of 
the  conducting  medium  the  second  term  in  the  bracket  is  small 
compared  to  the  constant  ‘3’  and  can  be  ignored.  Given  this,  the 
right  hand  side  of  this  equation  can  be  rewritten  as 


-MM  -Mo  p 


-sin  ^  cos  (40) 


£i  ..  VnM-e 

2  K  p ’ 

For  the  up-over-and-down  process,  ,  as  illustrated  in  Fig.  6, 
the  wave  induced  by  the  dipole  (HED)  first  propagates  upward 
(Part  1)  and  crosses  the  interface  into  the  free  space  region.  It  is 
attenuated  by  a  transmission  coefficient  Tup.  Second,  the  wave 
spreads  out  horizontally  (Part  II)  along  the  interface.  Note  that 
since  the  upper  medium  is  free  space  and  for  the  whole  range  of 


p  (100  -  10000m)  considered  here,  |/^p|«I,  the  fields  are 
quasi-static.  Finally,  at  the  position  on  the  interface  right  above 
P,  the  wave  crosses  the  interface  again  after  attenuation  by  a 
transmission  coefficient  T<jown  and  propagates  vertically  down 
(Part  III)  to  the  observation  point. 
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Fig.  6.  Illustration  oflhc  path  of  the  up-over-and-down  process 


To  identify  (40)  as  an  up-over-and-down  process,  it  is 
important  to  show  how  the  different  terms  have  functional 
dependencies  that  are  characteristic  to  different  portions  of  the 
process.  For  example,  it  will  be  shown  first  that  the  ‘up’  term  of 
(40)  represents  near  field  vertical  propagation  of  magnetic  field 
from  an  HED  in  a  homogeneous  medium  and  second,  why  this 
term  would  be  expected  to  be  a  magnetic  field  rather  than  an 
electric  field. 

To  begin,  it  is  helpful  to  understand  how  a  quasi-static  field 
in  the  air  can  be  excited  by  the  buried  HED.  This  can  be  done 
by  recognizing  that  the  HED  generates  tangential  electric  and 
magnetic  fields  just  above  the  air-water  interface  and  that  these 
fields  can  be  treated  as  equivalent  magnetic  and  electric  current 
sources  respectively  in  the  air  region  [7].  If,  further,  the  water  is 
replaced  by  a  perfect  magnetic  conductor  (PMC),  then  the  only 
remaining  source  is  the  electric  field  sources  that  originate  in 
the  incident  magnetic  field  and  are  doubled  in  value  by  the 
imaging  effect  [7].  The  use  of  image  theory  allows  the 
calculations  of  fields  to  be  done  in  a  homogeneous  air  medium 
These  fields  will  be  quasi-static  since  kop  «  1.  The  original 
case  for  HED  buried  in  PMC  is  shown  on  the  left  side  of  Fig.  7. 
The  tangential  magnetic  field  generated  by  the  HED  just  above 
the  air-water  interface  only  has  an  x  component,  H^nc  r.  The 

fields  in  the  air  region  can  be  obtained  by  the  equivalent  electric 
surface  current  source,  J  (A/m),  which  is  shown  on  the  right 


side  of  Fig.  7. 
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Fig.  7.  Equivalent  surface  current  source  for  the  HED  buried  in  PMC 
J  on  the  interface  (z  =  0)  can  be  determined  by 

7,=2n*Hl  =  2ayHl1  (41) 

where 
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/OO-'W 


Idl-he **' 
4 x  (r')J 


The  primed  variable  rl  represents  the  distance  from  the  dipole 
to  a  point  on  the  source  surface  and  r'  -  ^(p')2 +V  •  TupHx  is 

the  transmission  coefficient  for  jc  component  of  the  magnetic 
field. 

Tup.H*  (42) 

Integrating  the  source  current  over  the  entire  surface  gives 

a  =  f  f  J'P'tW  <43) 


where  Q  has  the  same  dimension  as  a  electric  dipole  moment 

(A-m).  Therefore,  it  can  be  treated  as  an  equivalent  dipole 
moment  to  replace  the  original  HED  if/?  is  much  larger  than  the 
size  of  the  source.  This  can  be  understood  as  follows.  The  HED 
is  first  replaced  by  the  surface  current  J  %  on  the  interface,  then 

the  distributed  surface  current  is  integrated  into  a  new 
horizontal  (y  orientation)  electric  dipole  just  above  the  interface. 
It  is  equivalent  to  say  that  the  original  HED  is  shifted  onto  the 
air-water  interface  with  a  change  in  magnitude.  The  new 
equivalent  HED  magnitude  can  be  determined  by  using  (41)  in 
the  integration  in  (43). 

Q,=ayIdle-*'XP„x  (44) 

The  exponential  term  in  (44)  is  obtained  by  assuming  that  k\r' 
varies  little  over  the  source  area  and 

e'Art  *e~Jk'h.  (45) 

As  shown  in  (4 1 ),  J s  is  dependent  on  the  depth  of  the  HED, 

h  It  is  proportional  to  Mh2  when  near  to  the  origin.  Beyond 
several  multiples  of  h  from  the  origin,  the  current  becomes  very 
small  and  negligible.  This  makes  the  approximation  in  (45) 
reasonable  because  within  several  multiples  of  h  the 
assumption  that  k\r'  varies  little  is  valid.  Fig.  8  shows  the 
magnitude  distribution  of  J  % . 
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Fig  8  Normalized  magnilude  distribution  of  the  surface  current  (all  the 
magnitudes  arc  normalized  by  the  maximum  of  the  curreni) 

The  magnitudes  in  the  figure  are  normalized  by  the 
magnitude  of  the  current  at  the  origin  (i.e.,  the  maximum  of  the 
current).  It  can  be  shown  that  there  is  an  effective  area  beyond 
which  the  effect  of  the  source  current  can  be  ignored  The 


effective  area  is  proportional  to  h 2.  Since  the  current  decays 
vertically  as  \fh2  and  the  effective  area  grows  with  h\  the 
integral  over  the  surface  will  be  independent  of  h.  Thus  the  new 
dipole  moment  as  shown  in  (44)  is  independent  of  h  except  for 
the  exponential  term.  The  process  shown  by  (41)  to  (45) 
describes  the  ‘up’  (Part  I  in  Fig.  6)  part  of  the  propagation. 

Given  the  new  equivalent  HED,  the  field  in  the  air  region  can 
be  determined  by  treating  the  new  HED  as  a  dipole  in  free  space. 
The  jc  component  of  the  electric  field  in  air  on  z  =  0  plane  is 

E°xJp,4>,  °)  =  ~  ~~  j— sin  <f> cos <t> 

4*o  P 

where  Qs  is  the  new  dipole  moment  and  the  subscript  V  stands 
for  the  field  obtained  by  the  up-over-and-down  interpretation. 
Replacing  Qs  by  (44)  results  in 

n  3  jrjJJl  ■  e-jk,he~jl"pTu 

El  ( p,  1 0)  =  -  - —— - affiL  Sin  ft  cos  <f>  (46) 

4*o  P 

This  is  a  field  propagating,  in  air,  over  the  air-water  interface 
and  the  ‘over’  (Part  II  in  Fig.  6)  part  of  the  propagation. 

The  field  shown  in  (46)  travels  from  the  interface  down  to 
the  observation  point  P  in  the  same  way  a  plane  wave  does. 

This  is  valid  because  \kcp\«l.  An  exponential  term,  C  , 
indicating  the  propagating  pattern  of  a  plane  wave  in  the  ‘-z’ 
direction  is  then  added  in  the  expression  of  the  field  at  P  to 
represent  the  ‘down’  (Part  111  in  Fig.  6)  part  of  the  propagation 
Finally,  the  expression  of  electric  field  component  Ex  at  the 
observation  point  is 

e\,Xp^  °)  = 


4*„ 


-sin^cos^T^.e 


V*l* 


(47) 


P 

which  shows  a  complete  process  of  the  up-over-and-down 
propagation.  The  transmission  coefficient  TdoNV11(px  in  (47)  is  for 
the  x  component  of  E  field  when  it  travels  down  and  crosses  the 
interface.  Since  there  is  no  reflected  wave  in  the  upper  region 
and  the  tangential  fields  arc  continuous  across  the  boundary, 

Tdown.Ex  —  1  • 

Comparing  (47)  to  (40),  there  is  a  constant  difference 
between  the  two  expressions  that 


E'xJp,t,z)  =  -2^TupMx  ■  E\(p,<f>,z) 

*o  n  i 

The  reason  why  the  difference  appears  is  still  not  identified  at 
this  point  But  if  Tup  jix  is  selected  as 

T  =2^A  =  3£°  (48) 

upJ  n0  *,  < 

then  the  simple  model  based  on  the  physical  interpretation  of 
the  propagation  process  matches  the  derived  expression. 


vill.  Conclusions 

A  method  to  simplify  the  Sommerfeld  integration  is 
proposed  in  this  paper.  Using  this  method  a  non-integral,  far 
field  approximation  of  the  E  &  H  field  due  to  a  horizontal 
electric  dipole  buried  in  a  conducting  half  space  is  obtained. 
Then  it  is  tested  that  this  approximation  works  very  well  when 
//,  z  « p,  |A/|  » |Atf|  and  \ktp\  »  I  and  the  error  is  less  than  1 0%. 
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This  approximation  reveals  the  up-over-and-down  behavior  of 
the  electromagnetic  wave  when  propagating  in  the  conducting 
half  space.  Finally  an  interpretation  based  on  the  surface 
equivalence  theorem  to  the  up-over-and-down  process  of  the 
wave  propagation  is  introduced. 
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Dual  Problem  Space  FDTD  Simulation  for 
Underwater  ELF  Applications 
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Abstract — The  finite-difTerenee  time-domain  (FDTD)  method  is 
being  used  to  simulate  extremely  low  frequencies  underwater.  In 
order  to  expand  the  potential  problem  space  without  reducing  the 
resolution  at  which  the  source  is  modeled,  a  near-to-far-field  trans¬ 
formation  method  has  been  developed. 

Index  Terms — Equivalent  sources,  finite-difference  time-domain 
(FDTD)  methods. 


I  Introduction 


Fig.  1  Two  problem  spaces  are  used  in  lhe  FDTD  simulation,  (a)  The  radiating 
source  is  modeled  in  the  near-held,  and  lhe  tangential  H-ficld  on  the  lhree-di- 
mcnsional  surface  are  calculated,  (b)  Using  lhe  tangeniial  H-fields  a  three-di 
mensional  source  is  impressed  in  the  far-field.  (a)  Near-field  problem  spate, 
(b)  Far-field  problem  space. 


UNDERWATER  mines  pose  the  greatest  threat  to  surface 
ships  [1].  These  mines  are  no  longer  restricted  to  direct 
contact,  but  can  detect  the  electromagnetic  (EM)  signature  of  a 
ship  and  launch  a  torpedo  from  hundreds  of  meters  away  [2], 
The  EM  radiation  of  a  surface  ship  tends  to  be  in  the  extremely 
low  frequency  (ELF)  range  and  can  propagate  underwater  over 
long  distances.  For  this  reason,  simulation  is  being  used  to  study 
underwater  EM  radiation. 

One  of  the  most  common  methods  used  in  EM  simulation  is 
the  finite-difference  time-domain  (FDTD)  method  [3],  [4].  It 
has  recently  been  shown  to  be  effective  at  ELF  frequencies  [5). 
However,  because  it  employs  a  uniform  grid,  one  is  always  left 
with  the  problem  of  choosing  a  cell  size  that  is  small  enough 
to  accurately  model  the  radiating  source,  but  large  enough 
to  model  an  extensive  far  field.  In  order  to  overcome  this,  a 
three-dimensional  near-to-far-field  transformation  has  been 
developed.  This  method  involves  two  separate  FDTD  simula¬ 
tions.  The  simulation  space  in  the  near-field  models  the  source, 
whether  it  is  a  ship’s  hull  or  an  antenna,  with  relatively  high 
resolution,  a  second  simulation  space  models  the  far-field  with 
larger  cells  in  order  to  model  propagation  hundreds  of  meters 
from  the  source.  The  transition  between  near-  and  far-field  is 
accomplished  by  applying  the  equivalence  principle  [6]. 

This  letter  is  arranged  as  follows.  The  use  of  the  equivalence 
principle  to  make  the  near-to-far-field  transition  is  described  in 
Section  II.  In  Section  III,  we  verify  the  accuracy  of  the  method 
by  comparison  with  an  analytic  method  based  on  Sommerfeld’s 
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half-space  method  [7]— [9].  In  Section  IV,  a  more  realistic  ex¬ 
ample  of  the  simulation  of  a  loop  antenna  in  a  lake  bed  is  pre¬ 
sented.  Section  V  summarizes  the  letter. 

II.  The  Near-To-Far-Field  Transformation 

Two  three-dimensional  FDTD  problem  spaces  are  utilized 
to  implement  the  near-to-far-field  transformation  (Fig.  1).  A 
problem  space  with  a  relatively  small  cell  size  (I  m3)  is  used 
to  model  the  source  [Fig.  1(a)].  Another  problem  space  with 
larger  cells  [13  m3)  is  used  to  model  the  far-field  [Fig.  1(b)]. 
The  ratio  of  13  to  1  between  far-  and  near-field  cell  sizes  re¬ 
sulted  in  near-  and  far-field  problem  spaces  of  about  the  same 
size,  which  seemed  to  be  the  optimum  case.  Each  problem  space 
is  surrounded  by  a  perfectly  matched  layer  (PML)  [10].  (A  new 
PML  for  ELF  frequencies  and  lossy  media  has  been  developed 
and  will  be  the  subject  of  a  future  paper.)  Each  problem  space 
contains  a  three-dimensional  transfer  surface  where  the  equiva¬ 
lence  principle  [6]  is  implemented.  On  a  surface,  the  source  is 
uniquely  specified  by  either  the  tangential  E-  or  H-fields.  We 
use  the  H-fields.  The  tangeniial  fields  calculated  on  the  transfer 
surface  in  the  near-field  are  impressed  on  the  transfer  surface  in 
the  far-field  to  form  the  far-field  source.  Since  the  far-field  cells 
are  13  times  larger  than  the  new  field  cells,  only  one  value  out 
of  13  in  the  near-field  is  needed  in  the  far-field.  This  method  is 
effective  even  when  the  medium  is  inhomogeneous,  as  will  be 
demonstrated  in  the  next  section.  There  is  no  coupling  from  the 
far-field  back  to  the  near-field.  This  is  illustrated  in  Fig.  2. 

Ill  Verification  of  the  Accuracy 

In  this  section,  the  results  of  the  near-to-far-field  transforma¬ 
tion  are  compared  to  analytic  results  calculated  using  Sommer¬ 
feld’s  half-space  (SHS)  problem.  SHS  problem  calculates  the 
resulting  fields  from  an  oscillating  dipole  near  a  plane  interface 
separating  two  homogeneous  half-space  regions,  as  illustrated 
in  Fig.  3.  This  method  is  well  described  in  the  literature  [7)-[9] 
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(a) 

Fig.  2.  (a)  Near-ficld  mesh  plot,  (b)  Far-field  mesh  plot. 
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Fig.  3  Diagram  of  the  three-dimensional  far-field  used  in  the  comparison  be¬ 
tween  the  FDTD  near-to-far-ficld  transformation  and  the  SHS  method.  Thex-di- 
rection  (not  shown)  is  1500  m.  The  source  is  a  magnetic  dipole  1  m  below  the 
air-water  interface  and  was  generated  in  the  near-field  (not  shown). 


TABLE  I 

Properties  op  the  Materials  Used  in  the  Simulations  described 
in  This  Letter  (5J 


Material  £r  O  (S  /  m) 

Air  1  0 

Lake  water  80  0  018 

Mud  40  0.002 

and  will  not  be  repeated  here.  In  Fig.  3,  the  upper  layer  is  air, 
the  middle  layer  is  water,  and  the  lower  level  is  mud.  The  water 
layer  in  the  middle  is  300  m  thick.  The  dielectric  properties  for 
water,  mud,  and  air  are  given  in  Table  I.  The  dipole  is  formed 
by  specifying  the  Hz-field  in  one  1-m3  cell  in  the  near-field 
problem  space.  The  monitor  lines  represent  the  places  where 
comparisons  between  the  methods  will  be  made.  Note  that  the 
FDTD  simulations  are  all  three-dimensional. 

Comparisons  at  200  and  1000  m  are  shown  in  Fig.  4(a)  and 
(b),  respectively.  The  simulation  required  30  OCX)  time-steps. 
The  amplitudes  in  each  figure  are  calculated  by  the  method  of 
two  equations,  two  unknowns  [11].  The  horizontal  coordinate  is 
the  distance  to  the  air/water  surface,  and  the  vertical  coordinate 
is  the  magnitude  of  the  field.  The  discrete  symbols  represent  the 
FDTD  calculations,  and  the  solid  lines  represent  the  calculations 
by  SHS  method.  Clearly,  the  results  of  the  comparisons  are  very 
good  in  all  cases. 


10  - * - * - * - - - * - 

0  50  100  150  200  250  300 
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Fig  4.  Comparisons  of  the  FDTD  simulations  (the  symbols)  and  the  SHS  cal¬ 
culations  (the  lines)  for  the  cell  size  ratio  of  13.  The  source  is  near  the  upper 
surface  of  the  water  layer.  The  comparisons  are  made  at  (a)  200  and  (b)  1000  m 
Z  represents  the  distance  from  the  surface,  (a)  Comparison  at  200  m.  (b)  Com¬ 
parison  at  1000  m. 


IV.  Example 

This  section  illustrates  the  use  of  FDTD  with  a  near-field, 
far-field  transformation  in  simulating  a  more  realistic  case.  One 
of  the  goals  of  the  simulation  is  to  verify  the  accuracy  of  the 


Authorized  licensed  use  limited  to  UNIVERSITY  OF  IDAHO.  Downloaded  on  Juty  27.2010  at  17:40:06  UTC  from  IEEE  Xplore  Restrictions  apply 


500 


IEEE  ANTENNAS  AND  WIRELESS  PROPAGATION  LETTERS.  VOL.  8,  2009 


4  meters 


h - 

< -  <- 


/ 


/ 


/ 


/ 


/ 


/ 


/ 


/  4  meters 


Ey  (source) 


(a) 


Hz(i+l/2j+l/2,k) 


/  i 

Hx(ij+l/2,k+l/2)  '  I 

/ 

/  i 

[  Ey(ij+l/2,k)j 

(7 


Hz(i-l/2j+l/2,k) 


Hx(ij+l/2,k-l/2) 


(b) 


Fig.  6  The  lake  bed  that  is  simulated  in  the  far-field.  The  cells  in  the  far  held 
arc  13  m3.  The  near-field  (the  enclosed  dashed  area)  contains  a  current  loop 
located  1  m  below  the  surface  of  the  water. 


Fig.  5.  The  simulation  of  a  4  m  x  4  m  current  loop  in  FDTD  as  generated  in 
the  near-field,  (a)  Setting  the  group  of  cells  at  the  respective  E-field  points  in  the 
grid  simulates  a  metal  loop,  (b)  Current  is  simulated  by  setting  an  Ev -field  to  a 
value  This  couples  to  the  surrounding  H-fields. 


modeling  with  measured  data  that  can  be  made  in  a  lake.  Be¬ 
cause  the  lake  bed  is  not  flat,  analytical  approaches  cannot  be 
used  to  solve  the  problem. 

The  antennas  that  will  be  used  in  the  transmission  are 
4  m  x  4  m  rectangular  current  loops.  Fig.  5(a)  illustrates  how 
this  is  simulated  in  the  XK-plane  in  the  FDTD  space.  Metal  can 
be  simulated  by  ensuring  that  an  E-field  is  zero  at  a  particular 
point  in  the  space.  Therefore,  using  cells  that  are  1  m3,  the 
metal  loop  antenna  is  simulated  at  the  corresponding  Ex  or  Ey 
positions,  as  shown  in  Fig.  5(a).  Since  the  radius  of  the  wire 
of  the  antennas  is  considerably  less  than  the  1-m  cell  size,  the 
thin  rod  approximation  [12]  is  used  to  model  the  wire  at  these 
positions 

In  FDTD,  a  current  cannot  be  simulated  directly,  but  it  can 
be  simulated  indirectly  by  using  Ampere's  circuit  law  [13]  and 
specifying  the  surrounding  H-fields 

I  =  <b—Hd\.  (1) 

J  Ho 

By  impressing  a  hard  source  on  one  of  the  Ey-fields,  a  value  is 
induced  on  the  surrounding  H-fields,  as  shown  in  Fig.  5(b).  This 
results  in  a  current  via  (1). 

A  model  of  the  lake  bed  is  created  for  the  far-field  domain 
(Fig.  6)  The  cells  are  13m3.  The  shape  of  the  lower  surface  of 
the  water  layer  shows  a  complex  geometry  structure  similar  to 
a  real  lake  bed. 


(a) 


(b) 

Fig.  7.  Results  of  the  simulation  illustrated  in  Fig.  6.  Z  represents  the  distance 
from  (he  surface,  (a)  The  results  at  200  m  (b)  The  results  at  1000  m. 

Fig.  7(a)  and(b)  shows  the  results  of  the  near-to-far-field  sim¬ 
ulation  using  the  current  loop  source  in  the  near-field  and  the 
lake  bed  in  the  far-field  at  200  and  1000  m  from  the  source  and 
at  three  different  frequencies,  10,  100,  and  1000  Hz.  These  sim¬ 
ulations  were  done  on  an  HP  DL140  GE  Quad  Core  and  re¬ 
quired  about  6  h.  Both  the  near-  and  far-field  problem  spaces 
were  120  cells  cubed. 

V.  Summary 

A  near-to-far-field  transformation  utilizing  the  equivalence 
principle  in  conjunction  with  the  FDTD  method  has  been  pre- 
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sented.  The  accuracy  of  this  method  was  confirmed  by  compar¬ 
isons  with  analytic  results  based  on  Sommerfeld’s  half-space 
method.  An  example  illustrating  the  flexibility  of  the  method  in 
simulating  a  realistic  problem  has  also  been  presented. 

The  method  presented  in  this  letter  substantially  extends  the 
range  of  FDTD  simulations  at  ELF  frequencies  for  the  purpose 
of  determining  the  vulnerability  of  surface  ships  to  electromag- 
netically  detonated  mines.  Accuracy  at  1  km  has  already  been 
confirmed,  and  it  is  hoped  that  the  development  of  a  new  under¬ 
water  PML  will  extend  this  range  to  3  km. 
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